56#include "CLHEP/Units/SystemOfUnits.h"
70#include <Eigen/StdVector>
81 return distsol.
first();
84 std::abs(distsol.
first()) < std::abs(distsol.
second()) ?
95 limitInversePValue(
double qOverP){
96 const double magnitude = std::abs(qOverP);
98 constexpr double maxP{100.*10e6*
MeV};
99 constexpr double minP{1.e-3*
MeV};
100 constexpr double lo {1./maxP};
101 constexpr double hi {1./minP};
102 const double limited = std::clamp(magnitude, lo, hi);
103 return std::copysign(limited, qOverP);
107 std::pair<const Trk::TrackParameters *, const Trk::TrackParameters *> getFirstLastIdPar(
const Trk::Track & track) {
113 while ((firstidpar ==
nullptr) && parit !=
track.trackParameters()->end()) {
115 ((**parit).covariance() !=
nullptr) &&
124 parit =
track.trackParameters()->end();
128 ((**parit).covariance() !=
nullptr) &&
133 }
while ((lastidpar ==
nullptr) && parit !=
track.trackParameters()->begin());
135 return std::make_pair(firstidpar, lastidpar);
150 std::abs(
a.parameters()[0] -
b.parameters()[0]) < e &&
151 std::abs(
a.parameters()[1] -
b.parameters()[1]) < e &&
152 std::abs(
a.parameters()[2] -
b.parameters()[2]) < e
163 calculateJac(Eigen::Matrix<double, 5, 5> &jac,
164 Eigen::Matrix<double, 5, 5> &out,
165 int jmin,
int jmax) {
169 out.block(0, 0, 4, jmin).setZero();
173 out.block(0, jmax + 1, 4, 5 - (jmax + 1)).setZero();
176 out(4, 4) = jac(4, 4);
186 std::pair<double, double> principalComponentAnalysis2x2(
const Amg::MatrixX & mat) {
187 const double trace =
mat(0, 0) +
mat(1, 1);
188 const double diagonalProduct =
mat(0, 0) *
mat(1, 1);
189 const double mat01Sq =
mat(0, 1) *
mat(0, 1);
190 const double discriminant = std::sqrt(trace * trace - 4. * (diagonalProduct - mat01Sq));
192 const double eigenValueSmall = 0.5 * (trace -
discriminant);
193 const double stereoAngle = 0.5 * std::asin(2 *
mat(0, 1) / (-discriminant));
195 return std::make_pair(eigenValueSmall, stereoAngle);
201 const std::string & t,
202 const std::string & n,
230 ATH_MSG_ERROR(
"Hole search requested but no boundary check tool provided.");
231 return StatusCode::FAILURE;
256 ATH_MSG_WARNING(
"FillDerivativeMatrix option selected, switching off acceleration!");
280 ATH_MSG_ERROR(
"Hole search requested but track summaries are disabled.");
281 return StatusCode::FAILURE;
286 return StatusCode::SUCCESS;
292 if (m_fit_status[
S_FITS] > 0) {
295 <<
" track fits failed because of a matrix inversion failure");
297 <<
" tracks were rejected by the outlier logic");
299 <<
" track fits failed because of a propagation failure");
301 <<
" track fits failed because of an invalid angle (theta/phi)");
303 <<
" track fits failed because the fit did not converge");
305 <<
" tracks did not pass the chi^2 cut");
307 <<
" tracks were killed by the energy loss update");
310 return StatusCode::SUCCESS;
315 std::unique_ptr<Track>
317 const EventContext& ctx,
323 ATH_MSG_DEBUG(
"--> entering GlobalChi2Fitter::fit(Track,Track,)");
331 !
cache.m_field_cache.solenoidOn() && !
cache.m_field_cache.toroidOn()
338 const Track *indettrack = firstismuon ? &intrk2 : &intrk1;
339 const Track *muontrack = firstismuon ? &intrk1 : &intrk2;
341 bool measphi =
false;
348 rot = &crot->rioOnTrack(0);
358 const double dotprod2 = measdir.dot(
361 if (std::abs(dotprod1) < .5 && std::abs(dotprod2) < .5) {
369 auto [firstidpar, lastidpar] = getFirstLastIdPar(*indettrack);
371 if ((firstidpar ==
nullptr) || (lastidpar ==
nullptr)) {
375 std::unique_ptr<const TrackParameters> parforcalo =
unique_clone(firstismuon ? firstidpar : lastidpar);
377 if (!
cache.m_field_cache.solenoidOn()) {
378 const AmgVector(5) & newpars = parforcalo->parameters();
380 parforcalo=parforcalo->associatedSurface().createUniqueTrackParameters(
381 newpars[0], newpars[1], newpars[2], newpars[3], 1 / 5000., std::nullopt
385 std::vector < MaterialEffectsOnTrack > calomeots;
388 calomeots =
m_calotool->extrapolationSurfacesAndEffects(
392 parforcalo->associatedSurface(),
405 if (calomeots.empty()) {
410 std::unique_ptr<Track> track;
413 cache.m_calomat =
false;
414 const bool tmp2 =
cache.m_extmat;
415 const bool tmp4 =
cache.m_idmat;
420 const double qoverpid = measperid !=
nullptr ? measperid->parameters()[
Trk::qOverP] : 0;
421 const double qoverpmuon = measpermuon !=
nullptr ? measpermuon->parameters()[
Trk::qOverP] : 0;
423 const AmgSymMatrix(5) * errmatid = measperid !=
nullptr ? measperid->covariance() :
nullptr;
424 const AmgSymMatrix(5) * errmatmuon = measpermuon !=
nullptr ? measpermuon->covariance() :
nullptr;
428 (errmatid !=
nullptr) &&
429 (errmatmuon !=
nullptr) &&
434 const double piderror = std::sqrt((*errmatid) (4, 4)) / (qoverpid * qoverpid);
435 const double pmuonerror = std::sqrt((*errmatmuon) (4, 4)) / (qoverpmuon * qoverpmuon);
436 const double energyerror = std::sqrt(
437 calomeots[1].energyLoss()->sigmaDeltaE() *
438 calomeots[1].energyLoss()->sigmaDeltaE() + piderror * piderror +
439 pmuonerror * pmuonerror
443 (std::abs(calomeots[1].energyLoss()->deltaE()) -
444 std::abs(1 / qoverpid) + std::abs(1 / qoverpmuon)
447 ATH_MSG_DEBUG(
"Changing from measured to parametrized energy loss");
452 parforcalo->associatedSurface(),
457 if (calomeots.empty()) {
465 bool firstfitwasattempted =
false;
468 if (!caloEntranceIsValid) {
473 (!
cache.m_field_cache.toroidOn() && !
cache.m_field_cache.solenoidOn()) ||
475 cache.m_getmaterialfromtrack &&
479 qoverpid * qoverpmuon > 0
484 if (
cache.m_fit_status[
S_FITS] == (
unsigned int) (nfits + 1)) {
485 firstfitwasattempted =
true;
490 (track ==
nullptr) &&
491 !firstfitwasattempted &&
492 (
cache.m_field_cache.toroidOn() ||
cache.m_field_cache.solenoidOn())
498 trajectory = trajectory2;
502 bool pseudoupdated =
false;
504 if (track !=
nullptr) {
505 for (std::unique_ptr<GXFTrackState> & pseudostate : trajectory.
trackStates()) {
506 if (pseudostate ==
nullptr) {
517 if ((pseudostate ==
nullptr) || pseudostate->fitQuality().chiSquared() < 10) {
522 const std::unique_ptr<const TrackParameters> updpar(
m_updator->removeFromState(
524 pseudostate->measurement()->localParameters(),
525 pseudostate->measurement()->localCovariance()
528 if (updpar ==
nullptr) {
533 covMatrix(0, 0) = 100;
535 std::unique_ptr<const PseudoMeasurementOnTrack> newpseudo = std::make_unique<const PseudoMeasurementOnTrack>(
539 std::move(covMatrix),
543 pseudostate->setMeasurement(std::move(newpseudo));
545 errors[0] = errors[2] = errors[3] = errors[4] = -1;
547 pseudostate->setMeasurementErrors(errors);
548 pseudoupdated =
true;
553 cache.m_matfilled =
true;
559 *track->perigeeParameters(),
564 cache.m_matfilled =
false;
570 if (track !=
nullptr) {
571 track->info().addPatternReco(intrk1.
info());
572 track->info().addPatternReco(intrk2.
info());
576 cache.m_calomat = tmp;
577 cache.m_extmat = tmp2;
578 cache.m_idmat = tmp4;
583 const EventContext& ctx,
585 const Track & intrk1,
586 const Track & intrk2,
588 std::vector<MaterialEffectsOnTrack> & calomeots
590 ATH_MSG_DEBUG(
"--> entering GlobalChi2Fitter::mainCombinationStrategy");
595 const Track *indettrack = firstismuon ? &intrk2 : &intrk1;
596 const Track *muontrack = firstismuon ? &intrk1 : &intrk2;
598 auto [tmpfirstidpar, tmplastidpar] = getFirstLastIdPar(*indettrack);
599 std::unique_ptr<const TrackParameters> firstidpar =
unique_clone(tmpfirstidpar);
600 std::unique_ptr<const TrackParameters> lastidpar =
unique_clone(tmplastidpar);
602 if ((firstidpar ==
nullptr) || (lastidpar ==
nullptr)) {
616 std::unique_ptr<const TrackParameters> tp_closestmuon =
nullptr;
618 while (closestmuonmeas ==
nullptr) {
619 closestmuonmeas =
nullptr;
622 if ((**tsosit).measurementOnTrack() !=
nullptr) {
623 closestmuonmeas = (**tsosit).measurementOnTrack();
625 if (thispar !=
nullptr) {
626 const AmgVector(5) & parvec = thispar->parameters();
628 parvec[0], parvec[1], parvec[2], parvec[3], parvec[4], std::nullopt
642 std::unique_ptr<const TrackParameters> tmppar;
645 if ((tp_closestmuon !=
nullptr) && msEntranceIsValid) {
650 std::unique_ptr<const std::vector<const TrackStateOnSurface *>> matvec;
652 if (tmppar !=
nullptr) {
653 const Surface & associatedSurface = tmppar->associatedSurface();
654 std::unique_ptr<Surface> muonsurf =
nullptr;
660 const double radius = cylbounds->
r();
662 muonsurf = std::make_unique<CylinderSurface>(trans, radius + 1, hlength);
666 const double newz = (
667 associatedSurface.
center().
z() > 0 ?
668 associatedSurface.
center().z() + 1 :
669 associatedSurface.
center().z() - 1
673 associatedSurface.
center().x(),
674 associatedSurface.
center().y(),
678 trans.translation() << newpos;
681 const double rmin = discbounds->
rMin();
682 const double rmax = discbounds->
rMax();
683 muonsurf = std::make_unique<DiscSurface>(trans, rmin, rmax);
687 if (muonsurf !=
nullptr) {
699 std::vector<const TrackStateOnSurface *> tmp_matvec;
701 if ((matvec !=
nullptr) && !matvec->empty()) {
702 tmp_matvec = *matvec;
703 delete tmp_matvec.back();
704 tmp_matvec.pop_back();
706 for (
auto & i : tmp_matvec) {
725 if (tmppar ==
nullptr) {
739 if (tmppar ==
nullptr) {
743 AmgVector(5) newpars = tmppar->parameters();
748 const double oldp = std::abs(1 / newpars[
Trk::qOverP]);
749 const double newp2 = oldp * oldp + (!firstismuon ? 2 : -2) * de * std::sqrt(mass * mass + oldp * oldp) + de * de;
756 tp_closestmuon=tmppar->associatedSurface().createUniqueTrackParameters(
757 newpars[0], newpars[1], newpars[2], newpars[3], newpars[4], std::nullopt
773 for (; itStates != endState; ++itStates) {
778 const bool tmpgetmat =
cache.m_getmaterialfromtrack;
780 if ((*itStates)->materialEffectsOnTrack() !=
nullptr) {
782 cache.m_extmat =
false;
784 cache.m_idmat =
false;
787 const auto *
const pBaseMEOT = (*itStates)->materialEffectsOnTrack();
791 const auto *
const pMEOT =
static_cast<const MaterialEffectsOnTrack *
>((*itStates)->materialEffectsOnTrack());
792 if ((pMEOT->scatteringAngles() ==
nullptr) or (pMEOT->energyLoss() ==
nullptr)) {
793 cache.m_getmaterialfromtrack =
true;
799 cache.m_getmaterialfromtrack = tmpgetmat;
806 trajectory.
trackStates().back()->setTrackParameters(
nullptr);
809 std::unique_ptr<const TrackParameters> firstscatpar;
810 std::unique_ptr<const TrackParameters> lastscatpar;
813 double newqoverpid = 0;
816 const double de = std::abs(calomeots[1].energyLoss()->deltaE());
817 const double sigmade = std::abs(calomeots[1].energyLoss()->sigmaDeltaE());
819 const double pbefore = std::abs(1 / firstidpar->parameters()[
Trk::qOverP]);
820 const double pafter = std::abs(1 / tp_closestmuon->parameters()[
Trk::qOverP]);
821 const double elosspull = (pbefore - pafter - de) / sigmade;
823 if (std::abs(elosspull) > 10) {
824 if (elosspull > 10) {
825 newqoverpid = 1 / (de + pafter + 10 * sigmade);
827 newqoverpid = 1 / (de + pafter - 10 * sigmade);
830 if (tp_closestmuon->parameters()[
Trk::qOverP] * newqoverpid < 0) {
834 const AmgVector(5) & newpar = firstidpar->parameters();
835 firstidpar=firstidpar->associatedSurface().createUniqueTrackParameters(
836 newpar[0], newpar[1], newpar[2], newpar[3], newqoverpid, std::nullopt
844 if (lastidpar ==
nullptr) {
850 *(firstismuon ? tp_closestmuon.get() : lastidpar.get()),
851 calomeots[0].associatedSurface(),
858 if (firstscatpar ==
nullptr) {
864 *(firstismuon ? firstidpar : tp_closestmuon),
865 calomeots[2].associatedSurface(),
872 if (lastscatpar ==
nullptr) {
876 std::optional<TransportJacobian> jac1;
877 std::optional<TransportJacobian> jac2;
878 std::unique_ptr<const TrackParameters> elosspar;
880 double firstscatphi = 0;
881 double secondscatphi = 0;
882 double firstscattheta = 0;
883 double secondscattheta = 0;
884 double muonscatphi = 0;
885 double muonscattheta = 0;
887 const TrackParameters *idscatpar = !firstismuon ? firstscatpar.get() : lastscatpar.get();
888 const TrackParameters *muonscatpar = firstismuon ? firstscatpar.get() : lastscatpar.get();
890 newqoverpid = idscatpar->parameters()[
Trk::qOverP];
892 const Amg::Vector3D calosegment = lastscatpar->position() - firstscatpar->position();
895 muonscattheta = calosegment.theta() - muonscatpar->parameters()[
Trk::theta];
898 for (
int i = 0; i < 2; i++) {
899 std::unique_ptr<const TrackParameters> tmpelosspar;
900 AmgVector(5) params1 = muonscatpar->parameters();
909 params1[0], params1[1], params1[2], params1[3], params1[4], std::nullopt
930 calomeots[1].associatedSurface(),
936 if ((tmpelosspar ==
nullptr) || (jac1 == std::nullopt)) {
945 const AmgVector(5) & newpars = tmpelosspar->parameters();
946 const std::unique_ptr<const TrackParameters> elosspar2(tmpelosspar->associatedSurface().createUniqueTrackParameters(
947 newpars[0], newpars[1], newpars[2], newpars[3], newqoverpid, std::nullopt
951 elosspar = std::move(tmpelosspar);
954 std::unique_ptr<const TrackParameters> scat2(
m_propagator->propagateParameters(
958 calomeots[0].associatedSurface() :
959 calomeots[2].associatedSurface(),
972 calomeots[0].associatedSurface() :
973 calomeots[2].associatedSurface(),
981 if ((scat2 ==
nullptr) || (jac2 == std::nullopt)) {
986 for (
int j = 0; j < 5; j++) {
987 for (
int k = 0; k < 5; k++) {
989 for (
int l = 0; l < 5; l++) {
990 jac3[j][k] += (*jac2) (j, l) * (*jac1) (l, k);
999 jac4(0, 0) = jac3[0][2];
1000 jac4(1, 1) = jac3[1][3];
1001 jac4(0, 1) = jac3[0][3];
1002 jac4(1, 0) = jac3[1][2];
1004 jac4 = jac4.inverse();
1016 discsurf =
static_cast<const Trk::DiscSurface *
>(&scat2->associatedSurface());
1018 if (cylsurf !=
nullptr) {
1019 dloc1 = -std::remainder(-dloc1, 2 *
M_PI * cylsurf->
bounds().
r());
1022 if (discsurf !=
nullptr) {
1023 dloc2 = -std::remainder(-dloc2, 2 *
M_PI);
1026 double dphi = jac4(0, 0) * dloc1 + jac4(0, 1) * dloc2;
1027 double dtheta = jac4(1, 0) * dloc1 + jac4(1, 1) * dloc2;
1033 muonscatphi += dphi;
1034 muonscattheta += dtheta;
1037 const double idscattheta = idscatpar->parameters()[
Trk::theta] - (scat2->parameters()[
Trk::theta] + dtheta);
1040 firstscatphi = muonscatphi;
1041 secondscatphi = idscatphi;
1042 firstscattheta = muonscattheta;
1043 secondscattheta = idscattheta;
1045 firstscatphi = -idscatphi;
1046 secondscatphi = -muonscatphi;
1047 firstscattheta = -idscattheta;
1048 secondscattheta = -muonscattheta;
1051 if (i == 1 &&
cache.m_field_cache.toroidOn() && !firstismuon) {
1052 AmgVector(5) params2 = scat2->parameters();
1060 firstscatpar=scat2->associatedSurface().createUniqueTrackParameters(
1061 params2[0], params2[1], params2[2], params2[3], params2[4], std::nullopt
1063 idscatpar = firstscatpar.get();
1067 *
cache.m_caloEntrance,
1071 if (startPar !=
nullptr) {
1080 rot.col(2) = trackdir;
1083 trans.linear().matrix() << rot;
1084 trans.translation() << startPar->position() - .1 * trackdir;
1095 if (curvlinpar !=
nullptr) {
1096 startPar.reset(curvlinpar);
1100 firstscatpar = std::move(scat2);
1104 std::unique_ptr<GXFMaterialEffects> firstscatmeff = std::make_unique<GXFMaterialEffects>(calomeots[0]);
1105 std::unique_ptr<GXFMaterialEffects> elossmeff = std::make_unique<GXFMaterialEffects>(calomeots[1]);
1106 std::unique_ptr<GXFMaterialEffects> secondscatmeff = std::make_unique<GXFMaterialEffects>(calomeots[2]);
1108 const double pull1 = std::abs(firstscatphi / firstscatmeff->sigmaDeltaPhi());
1109 const double pull2 = std::abs(secondscatphi / secondscatmeff->sigmaDeltaPhi());
1112 for (
auto & i : tmp_matvec) {
1117 firstscatmeff->setScatteringAngles(firstscatphi, firstscattheta);
1118 secondscatmeff->setScatteringAngles(secondscatphi, secondscattheta);
1121 elossmeff->setdelta_p(1000 * (lastscatpar->parameters()[
Trk::qOverP] - newqoverpid));
1123 elossmeff->setdelta_p(1000 * (newqoverpid - firstscatpar->parameters()[
Trk::qOverP]));
1126 elossmeff->setSigmaDeltaE(calomeots[1].energyLoss()->sigmaDeltaE());
1128 trajectory.
addMaterialState(std::make_unique<GXFTrackState>(std::move(firstscatmeff), std::move(firstscatpar)), -1);
1129 trajectory.
addMaterialState(std::make_unique<GXFTrackState>(std::move(elossmeff), std::move(elosspar)), -1);
1130 trajectory.
addMaterialState(std::make_unique<GXFTrackState>(std::move(secondscatmeff), std::move(lastscatpar)), -1);
1133 for (
auto & i : tmp_matvec) {
1140 if (startPar ==
nullptr) {
1145 (pull1 > 5 || pull2 > 5) &&
1151 bool largegap =
false;
1152 double previousz = 0;
1154 for (itStates2 = beginStates2; itStates2 != endState2; ++itStates2) {
1159 if (meff !=
nullptr) {
1172 cache.m_extmat =
false;
1174 cache.m_idmat =
false;
1180 !(itStates2 == beginStates2 || itStates2 == beginStates2 + 1) &&
1189 itStates2 == endState2 - 1 &&
1192 std::abs(tpar->
position().z()) < 13000
1194 std::unique_ptr<const TrackParameters> pseudopar(tpar->
clone());
1196 covMatrix(0, 0) = 100;
1198 std::unique_ptr<const PseudoMeasurementOnTrack> newpseudo = std::make_unique<const PseudoMeasurementOnTrack>(
1200 std::move(covMatrix),
1201 pseudopar->associatedSurface()
1204 std::unique_ptr<GXFTrackState> pseudostate = std::make_unique<GXFTrackState>(std::move(newpseudo), std::move(pseudopar));
1208 errors[0] = errors[2] = errors[3] = errors[4] = -1;
1211 pseudostate->setMeasurementErrors(errors);
1218 std::abs(trajectory.
trackStates().back()->position().z()) > 20000 &&
1219 std::abs(previousz) < 12000
1225 previousz = trajectory.
trackStates().back()->position().z();
1243 const EventContext& ctx,
1245 const Track & intrk1,
1246 const Track & intrk2,
1248 std::vector<MaterialEffectsOnTrack> & calomeots
1250 ATH_MSG_DEBUG(
"--> entering GlobalChi2Fitter::backupCombinationStrategy");
1253 const Track *indettrack = firstismuon ? &intrk2 : &intrk1;
1266 if (pParametersVector->size() > 1)
1267 firstidpar = (*pParametersVector)[1];
1269 firstidpar = pParametersVector->back();
1271 std::unique_ptr<const TrackParameters> lastidpar =
nullptr;
1272 if ((firstidpar !=
nullptr) && (
cache.m_caloEntrance !=
nullptr))
1276 if (lastidpar ==
nullptr) {
1277 lastidpar.reset(pParametersVector->back()->clone());
1280 std::unique_ptr < const TrackParameters > firstscatpar;
1281 std::unique_ptr < const TrackParameters > lastscatpar;
1282 std::unique_ptr < const TrackParameters > elosspar;
1287 lastidpar->position(),
1288 lastidpar->momentum(),
1290 lastidpar->position()
1297 calomeots[0].associatedSurface(),
1303 if (!firstscatpar) {
1307 const std::unique_ptr<const TrackParameters> tmppar(
1311 calomeots[1].associatedSurface(),
1324 const double oldp = std::abs(1 / tmppar->parameters()[
Trk::qOverP]);
1325 const double de = std::abs(calomeots[1].energyLoss()->deltaE());
1327 double newp2 = oldp * oldp - 2 * de * std::sqrt(mass * mass + oldp * oldp) + de * de;
1333 const double newqoverp =
sign / std::sqrt(newp2);
1335 const AmgVector(5) & pars = tmppar->parameters();
1338 tmppar->associatedSurface().createUniqueTrackParameters(
1339 pars[0], pars[1], pars[2], pars[3], newqoverp, std::nullopt
1345 calomeots[2].associatedSurface(),
1358 calomeots[2].associatedSurface(),
1372 calomeots[1].associatedSurface(),
1383 const double sign = (elosspar->parameters()[
Trk::qOverP] < 0) ? -1 : 1;
1384 const double newqoverp =
sign /
1385 (1. / std::abs(elosspar->parameters()[
Trk::qOverP]) +
1386 std::abs(calomeots[1].energyLoss()->deltaE()));
1388 const AmgVector(5) & pars = elosspar->parameters();
1390 std::unique_ptr<const TrackParameters>
const tmppar(
1391 elosspar->associatedSurface().createUniqueTrackParameters(
1392 pars[0], pars[1], pars[2], pars[3], newqoverp, std::nullopt
1399 calomeots[0].associatedSurface(),
1406 if (!firstscatpar) {
1411 for (; itStates != endState; ++itStates) {
1419 if ((*itStates)->materialEffectsOnTrack() !=
nullptr) {
1421 cache.m_idmat =
false;
1432 std::unique_ptr<GXFMaterialEffects> firstscatmeff = std::make_unique<GXFMaterialEffects>(calomeots[0]);
1433 std::unique_ptr<GXFMaterialEffects> elossmeff = std::make_unique<GXFMaterialEffects>(calomeots[1]);
1434 std::unique_ptr<GXFMaterialEffects> secondscatmeff = std::make_unique<GXFMaterialEffects>(calomeots[2]);
1438 sigmadp = calomeots[1].energyLoss()->sigmaDeltaE();
1439 elossmeff->setSigmaDeltaE(sigmadp);
1442 elossmeff->setdelta_p(dp);
1444 trajectory.
addMaterialState(std::make_unique<GXFTrackState>(std::move(firstscatmeff), std::move(firstscatpar)), -1);
1445 trajectory.
addMaterialState(std::make_unique<GXFTrackState>(std::move(elossmeff), std::move(elosspar)), -1);
1446 trajectory.
addMaterialState(std::make_unique<GXFTrackState>(std::move(secondscatmeff), std::move(lastscatpar)), -1);
1449 const Surface *triggersurf1 =
nullptr;
1450 const Surface *triggersurf2 =
nullptr;
1454 bool seenmdt =
false;
1455 bool mdtbetweenphihits =
false;
1459 itStates2 = (!firstismuon ? beginStates2 : endState - 1);
1460 itStates2 != (!firstismuon ? endState2 : beginStates - 1);
1461 (!firstismuon ? ++itStates2 : --itStates2)
1464 ((*itStates2)->measurementOnTrack() ==
nullptr) ||
1469 const auto *
const pMeasurement = (*itStates2)->measurementOnTrack();
1470 const Surface *surf = &pMeasurement->associatedSurface();
1474 if (isCompetingRIOsOnTrack) {
1476 rot = &crot->rioOnTrack(0);
1479 rot =
static_cast<const RIO_OnTrack *
>(pMeasurement);
1482 if ((rot !=
nullptr) &&
m_DetID->is_mdt(rot->
identify()) && (triggersurf1 !=
nullptr)) {
1486 (rot !=
nullptr) && (
1493 const double dotprod1 = measdir.dot(
Amg::Vector3D(0, 0, 1));
1496 const bool measphi = std::abs(dotprod1) <= .5 && std::abs(dotprod2) <= .5;
1500 (*itStates2)->trackParameters() !=
nullptr ?
1501 (*itStates2)->trackParameters()->position() :
1503 if (triggersurf1 !=
nullptr) {
1504 triggerpos2 = thispos;
1505 triggersurf2 = surf;
1507 mdtbetweenphihits =
true;
1510 triggerpos1 = thispos;
1511 triggersurf1 = surf;
1517 double mdttrig1 = 999999;
1518 double mdttrig2 = 999999;
1519 const Surface *mdtsurf1 =
nullptr;
1520 const Surface *mdtsurf2 =
nullptr;
1523 itStates2 = (!firstismuon ? beginStates2 : endState - 1);
1524 itStates2 != (!firstismuon ? endState2 : beginStates - 1);
1525 (!firstismuon ? ++itStates2 : --itStates2)
1527 const Surface *surf =
nullptr;
1529 ((*itStates2)->measurementOnTrack() !=
nullptr) &&
1532 surf = &(*itStates2)->measurementOnTrack()->associatedSurface();
1535 if (surf ==
nullptr) {
1538 const auto *
const pThisMeasurement = (*itStates2)->measurementOnTrack();
1542 if (isCompetingRioOnTrack) {
1544 rot = &crot->rioOnTrack(0);
1547 rot =
static_cast<const RIO_OnTrack *
>(pThisMeasurement);
1553 (*itStates2)->trackParameters() !=
nullptr ?
1554 (*itStates2)->trackParameters()->position() :
1555 pThisMeasurement->globalPosition();
1556 if (triggerpos1.mag() > 1 && (globpos - triggerpos1).mag() < mdttrig1) {
1557 mdttrig1 = (globpos - triggerpos1).
mag();
1560 if (triggerpos2.mag() > 1 && (globpos - triggerpos2).mag() < mdttrig2) {
1561 mdttrig2 = (globpos - triggerpos2).
mag();
1568 std::vector<GXFTrackState *> outlierstates;
1569 std::vector<GXFTrackState *> outlierstates2;
1571 outlierstates.reserve(10);
1572 outlierstates2.reserve(10);
1574 std::unique_ptr<PseudoMeasurementOnTrack> newpseudo;
1576 for (itStates2 = beginStates2; itStates2 != endState2; ++itStates2) {
1577 const auto *
const pMeasurement{(*itStates2)->measurementOnTrack()};
1579 const bool isStraightLine =
1580 pMeasurement !=
nullptr ?
1587 (newpseudo ==
nullptr) && (
1588 (itStates2 == beginStates2 || itStates2 == beginStates2 + 1) &&
1589 std::abs(pMeasurement->globalPosition().z()) < 10000
1592 std::unique_ptr<const TrackParameters> par2;
1593 if (((*itStates2)->trackParameters() !=
nullptr) && nphi > 99) {
1594 par2.reset((*itStates2)->trackParameters()->clone());
1606 if (par2 ==
nullptr) {
1610 covMatrix(0, 0) = 100;
1612 newpseudo = std::make_unique<PseudoMeasurementOnTrack>(
1614 std::move(covMatrix),
1615 par2->associatedSurface()
1618 std::unique_ptr<GXFTrackState> firstpseudo = std::make_unique<GXFTrackState>(std::move(newpseudo), std::move(par2));
1622 errors[0] = errors[2] = errors[3] = errors[4] = -1;
1625 firstpseudo->setMeasurementErrors(errors);
1626 firstpseudostate = firstpseudo.get();
1632 if (isPseudo && !firstismuon) {
1636 if ((**itStates2).materialEffectsOnTrack() !=
nullptr) {
1638 cache.m_idmat =
false;
1646 ((**itStates2).measurementOnTrack() !=
nullptr) &&
1647 &(**itStates2).measurementOnTrack()->associatedSurface() == triggersurf1 &&
1648 (mdtsurf1 !=
nullptr)
1650 std::unique_ptr<Amg::Transform3D> transf = std::make_unique<Amg::Transform3D>(mdtsurf1->
transform());
1652 transf->translation() << triggerpos1;
1655 covMatrix(0, 0) = 100;
1657 std::unique_ptr<const PseudoMeasurementOnTrack> newpseudo = std::make_unique<const PseudoMeasurementOnTrack>(
1661 std::unique_ptr<GXFTrackState> pseudostate1 = std::make_unique<GXFTrackState>(std::move(newpseudo),
nullptr);
1665 errors[0] = errors[2] = errors[3] = errors[4] = -1;
1668 pseudostate1->setMeasurementErrors(errors);
1669 outlierstates2.push_back(pseudostate1.get());
1674 ((**itStates2).measurementOnTrack() !=
nullptr) &&
1675 &(**itStates2).measurementOnTrack()->associatedSurface() == triggersurf2 &&
1676 mdtbetweenphihits &&
1677 (mdtsurf2 !=
nullptr)
1679 std::unique_ptr<Amg::Transform3D> transf = std::make_unique<Amg::Transform3D>(mdtsurf2->
transform());
1680 transf->translation() << triggerpos2;
1683 covMatrix(0, 0) = 100;
1685 std::unique_ptr<const PseudoMeasurementOnTrack> newpseudo = std::make_unique<const PseudoMeasurementOnTrack>(
1689 std::unique_ptr<GXFTrackState> pseudostate2 = std::make_unique<GXFTrackState>(std::move(newpseudo),
nullptr);
1693 errors[0] = errors[2] = errors[3] = errors[4] = -1;
1696 pseudostate2->setMeasurementErrors(errors);
1698 outlierstates2.push_back(pseudostate2.get());
1714 outlierstates.push_back(trajectory.
trackStates().back().get());
1722 Track *track =
nullptr;
1726 cache.m_matfilled =
true;
1727 const bool tmpacc =
cache.m_acceleration;
1728 cache.m_acceleration =
false;
1730 const std::unique_ptr<Trk::Track> tmp_track(
1732 cache.m_acceleration = tmpacc;
1734 cache.m_matfilled =
false;
1738 std::abs(trajectory.
residuals().tail<1>()(0) / trajectory.
errors().tail<1>()(0)) > 10
1744 if (firstpseudostate !=
nullptr) {
1747 covMatrix(0, 0) = 100;
1749 std::unique_ptr<const PseudoMeasurementOnTrack> newpseudo = std::make_unique<const PseudoMeasurementOnTrack>(
1751 std::move(covMatrix),
1759 for (
const auto & i : outlierstates2) {
1765 for (
const auto & i : outlierstates) {
1773 itStates = (firstismuon ? beginStates2 : endState - 1);
1774 itStates != (firstismuon ? endState2 : beginStates - 1);
1775 (firstismuon ? ++itStates : --itStates)
1788 cache.m_matfilled =
false;
1794 std::unique_ptr<Track>
1796 const EventContext& ctx,
1797 const Track& inputTrack,
1801 ATH_MSG_DEBUG(
"--> entering GlobalChi2Fitter::fit(Track,)");
1814 return std::unique_ptr<Track>(
1815 fitIm(ctx,
cache, inputTrack, runOutlier, matEffects));
1820 AlignmentCache& alignCache,
1821 const Track &inputTrack,
1828 alignCache.m_derivMatrix.reset();
1829 alignCache.m_fullCovarianceMatrix.reset();
1830 alignCache.m_iterationsOfLastFit = 0;
1833 fitIm(ctx,
cache, inputTrack, runOutlier, matEffects);
1834 if(newTrack !=
nullptr){
1835 if(
cache.m_derivmat.size() != 0)
1836 alignCache.m_derivMatrix = std::make_unique<Amg::MatrixX>(
cache.m_derivmat);
1837 if(
cache.m_fullcovmat.size() != 0)
1838 alignCache.m_fullCovarianceMatrix = std::make_unique<Amg::MatrixX>(
cache.m_fullcovmat);
1839 alignCache.m_iterationsOfLastFit =
cache.m_lastiter;
1846 const EventContext& ctx,
1848 const Track& inputTrack,
1853 ATH_MSG_DEBUG(
"--> entering GlobalChi2Fitter::fit(Track,,)");
1869 ATH_MSG_WARNING(
"Track without track parameters, cannot perform fit");
1877 if (minpar ==
nullptr) {
1881 const bool tmpgetmat =
cache.m_getmaterialfromtrack;
1887 cache.m_getmaterialfromtrack =
false;
1895 const Surface *firsthitsurf =
nullptr;
1896 const Surface *lasthitsurf =
nullptr;
1898 bool hasmuon =
false;
1899 bool iscombined =
false;
1900 bool seenphimeas =
false;
1904 for (; itStates != endState; ++itStates) {
1905 const auto *
const pMeasurement = (**itStates).measurementOnTrack();
1907 (pMeasurement ==
nullptr) &&
1908 ((**itStates).materialEffectsOnTrack() !=
nullptr) &&
1914 if (pMeasurement !=
nullptr) {
1915 const Surface *surf = &pMeasurement->associatedSurface();
1924 if (firsthitsurf ==
nullptr) {
1925 firsthitsurf = surf;
1928 if (
m_DetID->is_indet(hitid)) {
1933 if ((**itStates).trackParameters() !=
nullptr) {
1934 lastidpar = (**itStates).trackParameters();
1935 if (firstidpar ==
nullptr) {
1936 firstidpar = lastidpar;
1942 const double dotprod1 = measdir.dot(
Amg::Vector3D(0, 0, 1));
1944 if (std::abs(dotprod1) < .5 && std::abs(dotprod2) < .5) {
1946 if (std::abs(surf->
center().z()) > 13000) {
1949 if (surf->
center().perp() > 9000 && std::abs(surf->
center().z()) < 13000) {
1961 if (iscombined && seenphimeas && (phiem || phibo)) {
1972 (hasmuon ||
cache.m_acceleration)
1974 cache.m_matfilled =
true;
1977 if (firstidpar == lastidpar) {
1978 firstidpar = lastidpar =
nullptr;
1983 !
cache.m_matfilled &&
1988 (firstidpar !=
nullptr)
1997 const bool tmpacc =
cache.m_acceleration;
1999 const bool tmpsirecal =
cache.m_sirecal;
2000 std::unique_ptr<Track> tmptrack =
nullptr;
2004 cache.m_fiteloss =
true;
2005 cache.m_sirecal =
false;
2008 cache.m_asymeloss =
true;
2011 tmptrack.reset(
myfit(ctx,
cache, trajectory, *minpar,
false, matEffects));
2012 cache.m_sirecal = tmpsirecal;
2014 if (tmptrack ==
nullptr) {
2015 cache.m_matfilled =
false;
2016 cache.m_getmaterialfromtrack = tmpgetmat;
2017 cache.m_acceleration = tmpacc;
2018 cache.m_fiteloss = tmpfiteloss;
2023 bool isbrem =
false;
2025 unsigned int n_brem=0;
2027 for (std::unique_ptr<GXFTrackState> & state : trajectory.
trackStates()) {
2030 if (meff !=
nullptr) {
2036 const double p = 1. / std::abs(layerpars->parameters()[
Trk::qOverP] - .0005 * meff->
delta_p());
2038 std::optional<Amg::Vector2D> locpos(state->associatedSurface().globalToLocal(layerpars->
position()));
2039 const Amg::Vector3D layerNormal(state->associatedSurface().normal(*locpos));
2040 double costracksurf = 1.;
2042 costracksurf = std::abs(layerNormal.dot(layerpars->
momentum().unit()));
2044 const double oldde = meff->
deltaE();
2046 std::unique_ptr<EnergyLoss> eloss;
2047 double sigmascat = 0;
2049 if (matprop !=
nullptr) {
2050 eloss = std::make_unique<EnergyLoss>(
2051 m_elosstool->energyLoss(*matprop, p, 1. / costracksurf,
2053 sigmascat = std::sqrt(
m_scattool->sigmaSquare(*matprop, p, 1. / costracksurf, matEffects));
2055 if (eloss !=
nullptr) {
2060 sigmascat = std::sqrt(
m_scattool->sigmaSquare(tmpprop, p, 1. / costracksurf, matEffects));
2064 sigmascat / std::sin(layerpars->parameters()[
Trk::theta]),
2078 }
else if (eloss !=
nullptr) {
2089 refpars[0], refpars[1], refpars[2], refpars[3], refpars[4], std::nullopt
2093 cache.m_matfilled =
true;
2097 trajectory.
brems().clear();
2099 trajectory.
brems().resize(1);
2100 trajectory.
brems()[0] = bremdp;
2103 cache.m_asymeloss =
false;
2111 std::unique_ptr<Track> track(
myfit(ctx,
cache, trajectory, *minpar, runOutlier, matEffects));
2113 bool pseudoupdated =
false;
2115 if ((track !=
nullptr) && hasid && hasmuon) {
2116 for (std::unique_ptr<GXFTrackState> & pseudostate : trajectory.
trackStates()) {
2118 (pseudostate ==
nullptr) ||
2120 pseudostate->fitQuality().chiSquared() < 10
2126 const std::unique_ptr<const TrackParameters> updpar(
m_updator->removeFromState(
2128 pseudostate->measurement()->localParameters(),
2129 pseudostate->measurement()->localCovariance()
2132 if (updpar ==
nullptr) {
2137 covMatrix(0, 0) = 100;
2139 std::unique_ptr<const PseudoMeasurementOnTrack> newpseudo = std::make_unique<const PseudoMeasurementOnTrack>(
2141 std::move(covMatrix),
2145 pseudostate->setMeasurement(std::move(newpseudo));
2147 errors[0] = errors[2] = errors[3] = errors[4] = -1;
2149 pseudostate->setMeasurementErrors(errors);
2150 pseudoupdated =
true;
2153 if (pseudoupdated) {
2155 cache.m_matfilled =
true;
2156 track.reset(
myfit(ctx,
cache, trajectory, *track->perigeeParameters(),
false,
muon));
2157 cache.m_matfilled =
false;
2161 cache.m_matfilled =
false;
2162 cache.m_getmaterialfromtrack = tmpgetmat;
2163 cache.m_acceleration = tmpacc;
2164 cache.m_fiteloss = tmpfiteloss;
2166 if (track !=
nullptr) {
2169 track->info().addPatternReco(old_info);
2172 return track.release();
2175 std::unique_ptr<Track>
2177 const EventContext& ctx,
2183 ATH_MSG_DEBUG(
"--> entering GlobalChi2Fitter::fit(PRDS,TP,)");
2186 for (
const auto *prd : prds) {
2187 const Surface & prdsurf = (*prd).detectorElement()->surface((*prd).identify());
2192 plsurf =
static_cast < const PlaneSurface *
>(&prdsurf);
2199 if ((slsurf ==
nullptr) && (plsurf ==
nullptr)) {
2200 ATH_MSG_ERROR(
"Surface is neither PlaneSurface nor StraightLineSurface!");
2205 }
else if (slsurf !=
nullptr) {
2214 }
else if (plsurf !=
nullptr) {
2215 if (param.covariance() !=
nullptr) {
2237 if (rot !=
nullptr) {
2238 rots.push_back(rot);
2242 std::unique_ptr<Track> track =
2243 fit(ctx, rots, param, runOutlier, matEffects);
2245 for (
const auto *rot : rots) {
2252 std::unique_ptr<Track>
2254 const EventContext& ctx,
2255 const Track& inputTrack,
2260 ATH_MSG_DEBUG(
"--> entering GlobalChi2Fitter::fit(Track,Meas'BaseSet,,)");
2275 if (minpar ==
nullptr) {
2285 const bool old_reintoutl =
cache.m_reintoutl;
2286 cache.m_reintoutl =
false;
2287 const bool tmpasymeloss =
cache.m_asymeloss;
2290 cache.m_asymeloss =
true;
2293 for (; itStates != endState; ++itStates) {
2298 (trajectory.
trackStates().back()->materialEffects() !=
nullptr) &&
2299 trajectory.
trackStates().back()->materialEffects()->sigmaDeltaE() > 50.001
2301 trajectory.
trackStates().back()->materialEffects()->setKink(
true);
2305 cache.m_reintoutl = old_reintoutl;
2307 for (
const auto & measBase : addMeasColl) {
2308 if (measBase ==
nullptr) {
2309 ATH_MSG_WARNING(
"There is an empty MeasurementBase object in the track! Skip this object..");
2317 std::unique_ptr<Track> track(
myfit(ctx,
cache, trajectory, *minpar, runOutlier, matEffects));
2318 cache.m_asymeloss = tmpasymeloss;
2320 if (track !=
nullptr) {
2321 const double oldqual =
2326 const double newqual =
2327 track->fitQuality()->numberDoF() != 0 ?
2328 track->fitQuality()->chiSquared() / track->fitQuality()->numberDoF() :
2331 if (
m_extensioncuts && runOutlier && newqual > 2 && newqual > 2 * oldqual) {
2335 track.reset(
nullptr);
2339 if (track !=
nullptr) {
2348 std::unique_ptr<Track>
2350 const EventContext& ctx,
2356 ATH_MSG_DEBUG(
"--> entering GlobalChi2Fitter::fit(Track,PRDS,)");
2360 for (
const auto *prd : prds) {
2361 const Surface & prdsurf = (*prd).detectorElement()->surface((*prd).identify());
2364 std::unique_ptr<const TrackParameters>
const trackparForCorrect(
2380 rot =
m_ROTcreator->correct(*prd, *trackparForCorrect, ctx);
2383 if (rot !=
nullptr) {
2384 rots.push_back(rot);
2388 std::unique_ptr<Track> track =
fit(ctx,intrk, rots, runOutlier, matEffects);
2390 for (
const auto *rot : rots) {
2398 const EventContext& ctx,
2404 ATH_MSG_DEBUG(
"--> entering GlobalChi2Fitter::fit(Meas'BaseSet,,)");
2417 for (
const auto *itSet : rots) {
2418 if (itSet ==
nullptr) {
2419 ATH_MSG_WARNING(
"There is an empty MeasurementBase object in the track! Skip this object..");
2425 std::unique_ptr<const TrackParameters> startpar(param.
clone());
2428 matEffects ==
muon &&
2431 cache.m_matfilled =
true;
2434 myfit(ctx,
cache, trajectory, *startpar, runOutlier, matEffects);
2436 cache.m_matfilled =
false;
2450 covMatrix(0, 0) = 100;
2452 std::unique_ptr<const PseudoMeasurementOnTrack> newpseudo = std::make_unique<const PseudoMeasurementOnTrack>(
2454 std::move(covMatrix),
2458 trajectory.
trackStates().front()->setMeasurement(std::move(newpseudo));
2463 covMatrix(0, 0) = 100;
2465 std::unique_ptr<const PseudoMeasurementOnTrack> newpseudo = std::make_unique<const PseudoMeasurementOnTrack>(
2467 std::move(covMatrix),
2471 trajectory.
trackStates().back()->setMeasurement(std::move(newpseudo));
2478 refpars[0], refpars[1], refpars[2], refpars[3], refpars[4], std::nullopt
2483 myfit(ctx,
cache, trajectory, *startpar, runOutlier, matEffects);
2485 cache.m_matfilled =
true;
2494 refpars[0], refpars[1], refpars[2], refpars[3], refpars[4], std::nullopt
2501 firstpar = trajectory.
trackStates().front()->trackParameters();
2504 covMatrix(0, 0) = 100;
2506 std::unique_ptr<const PseudoMeasurementOnTrack> newpseudo = std::make_unique<const PseudoMeasurementOnTrack>(
2508 std::move(covMatrix),
2512 trajectory.
trackStates().front()->setMeasurement(std::move(newpseudo));
2514 errors[0] = errors[2] = errors[3] = errors[4] = -1;
2516 trajectory.
trackStates().front()->setMeasurementErrors(errors);
2520 lastpar = trajectory.
trackStates().back()->trackParameters();
2523 covMatrix(0, 0) = 100;
2525 std::unique_ptr<const PseudoMeasurementOnTrack> newpseudo = std::make_unique<const PseudoMeasurementOnTrack>(
2527 std::move(covMatrix),
2531 trajectory.
trackStates().back()->setMeasurement(std::move(newpseudo));
2533 errors[0] = errors[2] = errors[3] = errors[4] = -1;
2535 trajectory.
trackStates().back()->setMeasurementErrors(errors);
2539 std::unique_ptr<Track> track;
2541 if (startpar !=
nullptr) {
2542 track.reset(
myfit(ctx,
cache, trajectory, *startpar, runOutlier, matEffects));
2545 if (track !=
nullptr) {
2548 cache.m_matfilled =
false;
2565 ) &&
cache.m_getmaterialfromtrack
2575 std::unique_ptr<GXFMaterialEffects> newmeff;
2583 newmeff = std::make_unique<GXFMaterialEffects>(*meff);
2587 const double sigmascat = std::sqrt(
m_scattool->sigmaSquare(
2608 newmeff = std::make_unique<GXFMaterialEffects>(newmeot);
2624 ((**trajectory.
trackStates().rbegin()).trackParameters() !=
nullptr)
2626 const double delta_p = 1000 * (
2628 (**trajectory.
trackStates().rbegin()).trackParameters()->
2632 newmeff->setdelta_p(delta_p);
2668 seg =
static_cast<const Segment *
>(measbase);
2678 for (
int i = 0; i <
imax; i++) {
2681 std::unique_ptr<GXFTrackState> ptsos = std::make_unique<GXFTrackState>(
2682 std::unique_ptr<const MeasurementBase>(measbase2->
clone()),
2683 std::unique_ptr<const TrackParameters>(newtrackpar !=
nullptr ? newtrackpar->
clone() :
nullptr)
2686 double sinstereo = 0;
2688 errors[0] = errors[1] = errors[2] = errors[3] = errors[4] = -1;
2699 bool measphi =
false;
2702 bool rotated =
false;
2705 if (
m_DetID->is_pixel(hitid)) {
2707 }
else if (
m_DetID->is_sct(hitid)) {
2708 if (covmat.cols() != 1 && covmat(1, 0) != 0) {
2712 }
else if (
m_DetID->is_trt(hitid)) {
2722 }
else if (
m_DetID->is_mdt(hitid)) {
2724 }
else if (
m_DetID->is_tgc(hitid)) {
2729 }
else if (
m_DetID->is_csc(hitid)) {
2731 }
else if (
m_DetID->is_mm(hitid)) {
2733 }
else if (
m_DetID->is_stgc(hitid)) {
2739 const auto [covEigenValueSmall, covStereoAngle] = principalComponentAnalysis2x2(covmat);
2740 errors[0] = std::sqrt(covEigenValueSmall);
2741 sinstereo = std::sin(covStereoAngle);
2743 errors[0] = std::sqrt(covmat(0, 0));
2745 errors[1] = std::sqrt(covmat(1, 1));
2756 const double dotprod1 = measdir.dot(
Amg::Vector3D(0, 0, 1));
2758 if (std::abs(dotprod1) < .5 && std::abs(dotprod2) < .5) {
2766 int param_index = 0;
2768 errors[0] = std::sqrt(covmat(0, 0));
2773 errors[1] = std::sqrt(covmat(param_index, param_index));
2778 errors[2] = std::sqrt(covmat(param_index, param_index));
2783 errors[3] = std::sqrt(covmat(param_index, param_index));
2788 errors[4] = std::sqrt(covmat(param_index, param_index));
2809 ptsos->setMeasurementErrors(errors);
2810 ptsos->setSinStereo(sinstereo);
2811 ptsos->setMeasurementType(hittype);
2812 ptsos->setMeasuresPhi(measphi);
2814 if (isoutlier && !
cache.m_reintoutl) {
2833 if (tvol ==
nullptr) {
2837 const Trk::BinnedArray < Trk::Layer > *confinedLayers = tvol->
confinedLayers();
2840 if (confinedLayers !=
nullptr) {
2842 for (
const auto & layer : confinedLayers->
arrayObjects()) {
2844 if (layer !=
nullptr) {
2849 if ((layIndex.
value() == 0) || (layer->layerMaterialProperties() ==
nullptr)) {
2860 disclay =
static_cast<const DiscLayer *
>(layer);
2863 if (disclay !=
nullptr) {
2864 if (disclay->
center().z() < 0) {
2865 cache.m_negdiscs.push_back(disclay);
2867 cache.m_posdiscs.push_back(disclay);
2869 }
else if (cyllay !=
nullptr) {
2870 cache.m_barrelcylinders.push_back(cyllay);
2880 for (
size_t ib = 0 ; ib < bsurf.size(); ++ib) {
2881 const Layer *layer = bsurf[ib]->surfaceRepresentation().materialLayer();
2883 if (layer ==
nullptr)
continue;
2887 if ((layIndex.
value() == 0) || (layer->layerMaterialProperties() ==
nullptr)) {
2894 cylsurf =
static_cast<const CylinderSurface *
>(&layer->surfaceRepresentation());
2899 discsurf =
static_cast<const DiscSurface *
>(&layer->surfaceRepresentation());
2901 if (discsurf !=
nullptr) {
2903 discsurf->
center().z() < 0 &&
2904 std::find(
cache.m_negdiscs.begin(),
cache.m_negdiscs.end(), layer) ==
cache.m_negdiscs.end()
2906 cache.m_negdiscs.push_back(layer);
2908 discsurf->
center().z() > 0 &&
2909 std::find(
cache.m_posdiscs.begin(),
cache.m_posdiscs.end(), layer) ==
cache.m_posdiscs.end()
2911 cache.m_posdiscs.push_back(layer);
2914 (cylsurf !=
nullptr) &&
2915 std::find(
cache.m_barrelcylinders.begin(),
cache.m_barrelcylinders.end(), layer) ==
cache.m_barrelcylinders.end()
2917 cache.m_barrelcylinders.push_back(layer);
2920 if ((cylsurf ==
nullptr) && (discsurf ==
nullptr)) {
2927 if (confinedVolumes !=
nullptr) {
2928 for (
const auto & volume : confinedVolumes->
arrayObjects()) {
2929 if (volume !=
nullptr) {
2945 std::optional<std::pair<Amg::Vector3D, double>>
2958 const double * pos = parforextrap.
position().data();
2960 cache.m_field_cache.getFieldZR(pos, field);
2961 const double sinphi = std::sin(parforextrap.parameters()[
Trk::phi0]);
2962 const double cosphi = std::cos(parforextrap.parameters()[
Trk::phi0]);
2963 const double sintheta = std::sin(parforextrap.parameters()[
Trk::theta]);
2964 const double costheta = std::cos(parforextrap.parameters()[
Trk::theta]);
2968 const double r = (std::abs(currentqoverp) > 1e-10) ? -sintheta / (currentqoverp * 300. * field[2]) : 1e6;
2969 const double xc = parforextrap.
position().x() -
r * sinphi;
2970 const double yc = parforextrap.
position().y() +
r * cosphi;
2971 const double phi0 = std::atan2(parforextrap.
position().y() - yc, parforextrap.
position().x() - xc);
2972 const double z0 = parforextrap.
position().z();
2973 const double delta_s = (surf.
center().
z() -
z0) / costheta;
2974 const double delta_phi = delta_s * sintheta /
r;
2975 const double x = xc + std::abs(
r) * std::cos(
phi0 + delta_phi);
2976 const double y = yc + std::abs(
r) * std::sin(
phi0 + delta_phi);
2978 const double perp = intersect.perp();
2985 const double costracksurf = std::abs(costheta);
2987 return std::make_pair(intersect, costracksurf);
2990 std::optional<std::pair<Amg::Vector3D, double>>
3005 const double * pos = parforextrap.
position().data();
3007 cache.m_field_cache.getFieldZR(pos, field);
3008 const double sinphi = std::sin(parforextrap.parameters()[
Trk::phi0]);
3009 const double cosphi = std::cos(parforextrap.parameters()[
Trk::phi0]);
3010 const double sintheta = std::sin(parforextrap.parameters()[
Trk::theta]);
3011 const double costheta = std::cos(parforextrap.parameters()[
Trk::theta]);
3012 const double tantheta = std::tan(parforextrap.parameters()[
Trk::theta]);
3013 const double r = (std::abs(currentqoverp) > 1e-10) ? -sintheta / (currentqoverp * 300. * field[2]) : 1e6;
3014 const double xc = parforextrap.
position().x() -
r * sinphi;
3015 const double yc = parforextrap.
position().y() +
r * cosphi;
3016 const double phi0 = std::atan2(parforextrap.
position().y() - yc, parforextrap.
position().x() - xc);
3017 const double z0 = parforextrap.
position().z();
3018 const double d = xc * xc + yc * yc;
3019 const double rcyl = surf.
bounds().
r();
3020 double mysqrt = ((
r + rcyl) * (
r + rcyl) - d) * (d - (
r - rcyl) * (
r - rcyl));
3026 mysqrt = std::sqrt(mysqrt);
3027 double firstterm = xc / 2 + (xc * (rcyl * rcyl -
r *
r)) / (2 * d);
3028 double secondterm = (mysqrt * yc) / (2 * d);
3029 const double x1 = firstterm + secondterm;
3030 const double x2 = firstterm - secondterm;
3031 firstterm = yc / 2 + (yc * (rcyl * rcyl -
r *
r)) / (2 * d);
3032 secondterm = (mysqrt * xc) / (2 * d);
3033 const double y1 = firstterm - secondterm;
3034 const double y2 = firstterm + secondterm;
3037 const double dist1 = (
x - x1) * (
x - x1) + (
y - y1) * (
y - y1);
3038 const double dist2 = (
x - x2) * (
x - x2) + (
y - y2) * (
y - y2);
3040 if (dist1 < dist2) {
3048 const double phi1 = std::atan2(
y - yc,
x - xc);
3051 const double delta_z =
r * deltaphi / tantheta;
3052 const double z =
z0 + delta_z;
3063 const Amg::Vector3D trackdir(cos(phidir) * sintheta, std::sin(phidir) * sintheta, costheta);
3065 const double costracksurf = std::abs(normal.unit().dot(trackdir));
3067 return std::make_pair(intersect, costracksurf);
3074 std::vector<std::pair<const Layer *, const Layer *>> &layers,
3079 std::vector<std::unique_ptr<GXFTrackState>> & states = trajectory.
trackStates();
3080 std::vector<std::unique_ptr<GXFTrackState>> oldstates = std::move(states);
3083 states.reserve(oldstates.size() + layers.size());
3091 for (
int i = 0; i <= indexoffset; i++) {
3101 for (
int i = indexoffset + 1; i < (int) oldstates.size(); i++) {
3102 const double rmeas = oldstates[i]->position().perp();
3103 const double zmeas = oldstates[i]->position().z();
3111 while (layerindex < (
int) layers.size()) {
3113 double costracksurf = 0.0;
3114 const Layer *layer =
nullptr;
3123 if (layers[layerindex].first !=
nullptr) {
3127 layer = layers[layerindex].first;
3134 if (oldstates[i]->trackParameters() !=
nullptr) {
3135 const double rlayer = cylsurf->
bounds().
r();
3136 if (std::abs(rmeas - rlayer) < std::abs(parforextrap->
position().perp() - rlayer)) {
3137 parforextrap = oldstates[i]->trackParameters();
3147 std::tie(intersect, costracksurf) =
res.value();
3153 if (cylsurf->
bounds().
r() > rmeas)
break;
3154 }
else if (layers[layerindex].second !=
nullptr) {
3160 layer = layers[layerindex].second;
3163 if (oldstates[i]->trackParameters() !=
nullptr) {
3164 const double zlayer = discsurf->
center().z();
3165 if (std::abs(zmeas - zlayer) < std::abs(parforextrap->
position().z() - zlayer)) {
3166 parforextrap = oldstates[i]->trackParameters();
3171 std::tie(intersect, costracksurf) =
res.value();
3177 if (std::abs(discsurf->
center().z()) > std::abs(zmeas))
break;
3179 throw std::logic_error(
"Unhandled surface.");
3186 const MaterialProperties *matprop = layer->layerMaterialProperties()->fullMaterial(intersect);
3187 if (matprop ==
nullptr) {
3200 const double actualx0 = X0 / costracksurf;
3201 const double de = -std::abs(
3205 (
m_p != 0.0 ? std::abs(
m_p) : std::abs(1. / currentqoverp)),
3207 const double sintheta = std::sin(parforextrap->parameters()[
Trk::theta]);
3208 const double sigmascat = std::sqrt(
m_scattool->sigmaSquare(
3210 (
m_p != 0.0 ? std::abs(
m_p) : std::abs(1. / currentqoverp)),
3214 std::unique_ptr<GXFMaterialEffects> meff = std::make_unique<GXFMaterialEffects>();
3215 meff->setDeltaE(de);
3216 meff->setScatteringSigmas(sigmascat / sintheta, sigmascat);
3217 meff->setX0(actualx0);
3218 meff->setSurface(&layer->surfaceRepresentation());
3219 meff->setMaterialProperties(matprop);
3225 std::unique_ptr<EnergyLoss> eloss;
3228 eloss = std::make_unique<EnergyLoss>(
m_elosstool->energyLoss(
3230 (
m_p != 0.0 ? std::abs(
m_p) : std::abs(1. / currentqoverp)),
3235 if (eloss !=
nullptr) {
3236 meff->setSigmaDeltaE(eloss->sigmaDeltaE());
3241 meff->setDeltaE(-5);
3243 meff->setScatteringSigmas(0, 0);
3246 meff->setSigmaDeltaE(50);
3247 if (eloss !=
nullptr) {
3248 meff->setSigmaDeltaEPos(eloss->sigmaPlusDeltaE());
3249 meff->setSigmaDeltaENeg(eloss->sigmaMinusDeltaE());
3254 "X0: " << meff->x0() <<
" qoverp: " << currentqoverp <<
3255 " sigmascat " << meff->sigmaDeltaTheta() <<
" eloss: " << meff->deltaE() <<
3256 " sigma eloss: " << meff->sigmaDeltaE()
3263 std::unique_ptr<GXFTrackState> matstate = std::make_unique<GXFTrackState>(
3265 std::unique_ptr<const TrackParameters>()
3267 matstate->setPosition(intersect);
3283 std::vector<std::pair<const Layer *, const Layer *>> & layers,
3284 std::vector<std::pair<const Layer *, const Layer *>> & upstreamlayers,
3285 const std::vector<std::unique_ptr<GXFTrackState>> & oldstates,
3294 upstreamlayers.reserve(5);
3313 const double lastz = laststate->
position().z();
3314 const double lastr = laststate->
position().perp();
3320 const double tantheta = std::tan(refpar->parameters()[
Trk::theta]);
3321 const double slope = (tantheta != 0) ? 1 / tantheta : 0;
3327 std::vector < const Layer *>::const_iterator it;
3328 std::vector < const Layer *>::const_iterator itend;
3336 it =
cache.m_posdiscs.begin();
3337 itend =
cache.m_posdiscs.end();
3339 it =
cache.m_negdiscs.begin();
3340 itend =
cache.m_negdiscs.end();
3346 for (; it != itend; ++it) {
3351 if (std::abs((*it)->surfaceRepresentation().center().z()) > std::abs(lastz)) {
3359 const DiscBounds *discbounds =
static_cast<const DiscBounds *
> (&(*it)->surfaceRepresentation().bounds());
3364 if (discbounds->
rMax() < firstr || discbounds->
rMin() > lastr) {
3368 const double rintersect = firstr + ((*it)->surfaceRepresentation().center().
z() - firstz) / slope;
3371 rintersect < discbounds->rMin() - 50 ||
3372 rintersect > discbounds->
rMax() + 50
3382 if ((*it) == endlayer) {
3394 std::abs((*it)->surfaceRepresentation().center().z()) < std::abs(firstz) ||
3397 upstreamlayers.emplace_back((
Layer *)
nullptr, (*it));
3404 (*it) != startlayer &&
3405 (std::abs((*it)->surfaceRepresentation().center().z()) > std::abs(firstz2) ||
3406 (*it) == startlayer2)
3408 layers.emplace_back((
Layer *)
nullptr, (*it));
3416 for (
const auto *barrelcylinder :
cache.m_barrelcylinders) {
3420 if (barrelcylinder->surfaceRepresentation().bounds().r() > lastr) {
3427 const double zintersect = firstz + (barrelcylinder->surfaceRepresentation().bounds().
r() - firstr) * slope;
3429 if (std::abs(zintersect - barrelcylinder->surfaceRepresentation().center().z()) >
3430 ((
const CylinderSurface*)(&barrelcylinder->surfaceRepresentation()))->bounds().halflengthZ() + 50) {
3434 if (barrelcylinder == endlayer) {
3441 if (barrelcylinder->surfaceRepresentation().bounds().r() < firstr ||
3442 barrelcylinder == startlayer) {
3443 upstreamlayers.emplace_back(barrelcylinder, (
Layer*)
nullptr);
3446 if (barrelcylinder != startlayer &&
3447 (barrelcylinder->surfaceRepresentation().bounds().r() > firstr2 ||
3448 barrelcylinder == startlayer2)) {
3449 layers.emplace_back(barrelcylinder, (
Layer*)
nullptr);
3462 const EventContext& ctx,
3473 if (!caloEntranceIsValid) {
3482 cache.m_negdiscs.empty() &&
3483 cache.m_posdiscs.empty() &&
3484 cache.m_barrelcylinders.empty()
3499 cache.m_fastmat =
false;
3516 bool hasmat =
false;
3517 int indexoffset = 0;
3518 int lastmatindex = 0;
3519 std::vector<std::unique_ptr<GXFTrackState>> & oldstates = trajectory.
trackStates();
3533 for (
int i = 0; i < (int) oldstates.size(); i++) {
3534 if (oldstates[i]->materialEffects() !=
nullptr) {
3543 if (firstsistate ==
nullptr) {
3544 if (oldstates[i]->trackParameters() ==
nullptr) {
3545 std::unique_ptr<const TrackParameters> tmppar(
m_propagator->propagateParameters(
3548 oldstates[i]->associatedSurface(),
3555 if (tmppar ==
nullptr)
return;
3557 oldstates[i]->setTrackParameters(std::move(tmppar));
3559 firstsistate = oldstates[i].get();
3561 lastsistate = oldstates[i].get();
3569 if (lastsistate ==
nullptr) {
3570 throw std::logic_error(
"No track state");
3580 std::unique_ptr<const TrackParameters> tmppar(
m_propagator->propagateParameters(
3589 if (tmppar ==
nullptr)
return;
3601 indexoffset = lastmatindex;
3616 std::vector<std::pair<const Layer *, const Layer *>> layers;
3617 std::vector<std::pair<const Layer *, const Layer *>> & upstreamlayers = trajectory.
upstreamMaterialLayers();
3631 const EventContext& ctx,
3637 if (refpar2 ==
nullptr) {
3647 std::vector<std::unique_ptr<GXFTrackState>> & states = trajectory.
trackStates();
3648 std::vector<std::unique_ptr<GXFTrackState>> matstates;
3649 std::unique_ptr< const std::vector < const TrackStateOnSurface *>,
3650 void (*)(
const std::vector<const TrackStateOnSurface *> *) >
3652 bool matvec_used=
false;
3653 std::unique_ptr<TrackParameters> startmatpar1;
3654 std::unique_ptr<TrackParameters> startmatpar2;
3665 int npseudomuon1 = 0;
3666 int npseudomuon2 = 0;
3668 for (
auto & state : states) {
3674 if (firstidhit ==
nullptr) {
3683 if (firsthit ==
nullptr) {
3684 firsthit = state->measurement();
3685 if (
cache.m_acceleration) {
3686 if (tp ==
nullptr) {
3696 if (tp ==
nullptr) {
3700 state->setTrackParameters(std::unique_ptr<const TrackParameters>(tp));
3706 lasthit = state->measurement();
3712 if (firstidhit ==
nullptr) {
3713 firstidhit = state->measurement();
3716 if ((firstidpar ==
nullptr) && (tp !=
nullptr)) {
3720 lastidhit = state->measurement();
3721 if (tp !=
nullptr) {
3726 if (firstsiliconpar ==
nullptr) {
3727 firstsiliconpar = tp;
3729 lastsiliconpar = tp;
3741 if (firstmuonhit ==
nullptr) {
3742 firstmuonhit = state->measurement();
3743 if (tp !=
nullptr) {
3747 lastmuonhit = state->measurement();
3748 if (tp !=
nullptr) {
3754 if (meff->
deltaE() == 0) {
3755 if (firstcalopar ==
nullptr) {
3756 firstcalopar = state->trackParameters();
3758 lastcalopar = state->trackParameters();
3760 if (firstmatpar ==
nullptr) {
3761 firstmatpar = state->trackParameters();
3766 std::unique_ptr<TrackParameters> refpar;
3767 AmgVector(5) newpars = refpar2->parameters();
3774 newpars[0], newpars[1], newpars[2], newpars[3], newpars[4], std::nullopt
3777 if (firstmatpar !=
nullptr) {
3782 if ((startmatpar1 ==
nullptr) || ((firstidhit !=
nullptr) && (firstmuonhit !=
nullptr))) {
3786 const double mass = trajectory.
mass();
3787 if (mass > 200 *
MeV) {
3788 const AmgVector(5) & newpars = startmatpar2->parameters();
3789 const double oldp = std::abs(1 / newpars[
Trk::qOverP]);
3792 startmatpar2 = startmatpar2->associatedSurface().createUniqueTrackParameters(
3793 newpars[0], newpars[1], newpars[2], newpars[3],
3794 sign / std::sqrt(oldp * oldp + 2 * 100 *
MeV * std::sqrt(oldp * oldp + mass * mass) + 100 *
MeV * 100 *
MeV),
3799 AmgVector(5) newpars = startmatpar1->parameters();
3802 startmatpar1 = startmatpar1->associatedSurface().createUniqueTrackParameters(
3803 newpars[0], newpars[1], newpars[2], newpars[3], newpars[4], std::nullopt
3806 newpars = startmatpar2->parameters();
3809 startmatpar2 = startmatpar2->associatedSurface().createUniqueTrackParameters(
3810 newpars[0], newpars[1], newpars[2], newpars[3], newpars[4], std::nullopt
3818 refpar->momentum().unit()
3821 const double distance = getDistance(distsol);
3824 ATH_MSG_DEBUG(
"Obtaining upstream layers from Extrapolator");
3827 std::unique_ptr<const TrackParameters> tmppar;
3829 if (firstmuonhit !=
nullptr) {
3831 if (caloEntranceIsValid) {
3834 *
cache.m_caloEntrance,
3838 if (tmppar !=
nullptr) {
3839 destsurf = &tmppar->associatedSurface();
3844 if (matvec_used)
cache.m_matTempStore.push_back( std::move(matvec) );
3849 false, matEffects) );
3852 if (matvec && !matvec->empty()) {
3853 for (
int i = (
int)matvec->size() - 1; i > -1; i--) {
3858 std::unique_ptr<GXFMaterialEffects> meff = std::make_unique<GXFMaterialEffects>(*meot);
3859 const TrackParameters * newpars = (*matvec)[i]->trackParameters() !=
nullptr ? (*matvec)[i]->trackParameters()->clone() :
nullptr;
3860 meff->setSigmaDeltaE(0);
3861 matstates.push_back(std::make_unique<GXFTrackState>(
3863 std::unique_ptr<const TrackParameters>(newpars)
3876 refpar->momentum().unit()
3879 const double distance = getDistance(distsol);
3882 ATH_MSG_DEBUG(
"Obtaining downstream ID layers from Extrapolator");
3884 std::unique_ptr<const TrackParameters> tmppar;
3885 std::unique_ptr<Surface> calosurf;
3886 if (firstmuonhit !=
nullptr) {
3888 if (caloEntranceIsValid) {
3891 *
cache.m_caloEntrance,
3896 if (tmppar !=
nullptr) {
3900 cylcalosurf =
static_cast<const CylinderSurface *
>(&tmppar->associatedSurface());
3905 disccalosurf =
static_cast<const DiscSurface *
>(&tmppar->associatedSurface());
3907 if (cylcalosurf !=
nullptr) {
3910 const double radius = cylbounds.
r();
3912 calosurf = std::make_unique<CylinderSurface>(trans, radius - 1, hlength);
3913 }
else if (disccalosurf !=
nullptr) {
3914 const double newz = (
3915 disccalosurf->
center().
z() > 0 ?
3916 disccalosurf->
center().z() - 1 :
3917 disccalosurf->
center().z() + 1
3921 disccalosurf->
center().x(),
3922 disccalosurf->
center().y(),
3927 trans.translation() << newpos;
3930 const double rmin = discbounds->
rMin();
3931 const double rmax = discbounds->
rMax();
3932 calosurf = std::make_unique<DiscSurface>(trans, rmin, rmax);
3934 destsurf = calosurf.release();
3938 if (matvec_used)
cache.m_matTempStore.push_back( std::move(matvec) );
3940 ctx, *startmatpar2, *destsurf,
alongMomentum,
false, matEffects));
3941 matvec_used =
false;
3943 if (matvec && !matvec->empty()) {
3944 for (
const auto & i : *matvec) {
3950 std::unique_ptr<GXFMaterialEffects> meff = std::make_unique<GXFMaterialEffects>(*meot);
3956 meff->setDeltaE(-5);
3959 meff->setScatteringSigmas(0, 0);
3962 meff->setSigmaDeltaE(50);
3965 const TrackParameters * newparams = i->trackParameters() !=
nullptr ? i->trackParameters()->clone() :
nullptr;
3967 matstates.push_back(std::make_unique<GXFTrackState>(
3969 std::unique_ptr<const TrackParameters>(newparams)
3981 if (
cache.m_calomat && (firstmuonhit !=
nullptr) && (firstidhit !=
nullptr)) {
3984 std::vector<MaterialEffectsOnTrack> calomeots =
m_calotool->extrapolationSurfacesAndEffects(
3993 if (calomeots.empty()) {
3998 std::unique_ptr<const TrackParameters> prevtrackpars =
unique_clone(lastidpar);
3999 if (lasthit == lastmuonhit) {
4000 for (
int i = 0; i < (int) calomeots.size(); i++) {
4003 std::unique_ptr<const TrackParameters> layerpar(
m_propagator->propagateParameters(
4006 calomeots[i].associatedSurface(),
4013 if (layerpar ==
nullptr) {
4018 std::unique_ptr<GXFMaterialEffects> meff = std::make_unique<GXFMaterialEffects>(calomeots[i]);
4021 lastcalopar = layerpar.get();
4025 const double qoverp = layerpar->parameters()[
Trk::qOverP];
4026 double qoverpbrem = 0;
4030 (firstmuonpar !=
nullptr) &&
4031 std::abs(firstmuonpar->parameters()[
Trk::qOverP]) > 1.e-9
4033 qoverpbrem = firstmuonpar->parameters()[
Trk::qOverP];
4035 const double sign = (qoverp > 0) ? 1 : -1;
4036 qoverpbrem =
sign / (1 / std::abs(qoverp) - std::abs(calomeots[i].energyLoss()->deltaE()));
4039 const AmgVector(5) & newpar = layerpar->parameters();
4041 layerpar = layerpar->associatedSurface().createUniqueTrackParameters(
4042 newpar[0], newpar[1], newpar[2], newpar[3], qoverpbrem, std::nullopt
4044 meff->setdelta_p(1000 * (qoverpbrem - qoverp));
4047 matstates.push_back(std::make_unique<GXFTrackState>(
4049 std::unique_ptr<const TrackParameters>(layerpar !=
nullptr ? layerpar->clone() :
nullptr)
4051 prevtrackpars = std::move(layerpar);
4056 firsthit == firstmuonhit &&
4057 (!
cache.m_getmaterialfromtrack || lasthit == lastidhit)
4060 for (
int i = 0; i < (int) calomeots.size(); i++) {
4062 std::unique_ptr<const TrackParameters> layerpar(
m_propagator->propagateParameters(
4065 calomeots[i].associatedSurface(),
4072 if (layerpar ==
nullptr) {
4077 std::unique_ptr<GXFMaterialEffects> meff = std::make_unique<GXFMaterialEffects>(calomeots[i]);
4086 const double qoverpbrem = layerpar->parameters()[
Trk::qOverP];
4091 (lastmuonpar !=
nullptr) &&
4092 std::abs(lastmuonpar->parameters()[
Trk::qOverP]) > 1.e-9
4096 const double sign = (qoverpbrem > 0) ? 1 : -1;
4097 qoverp =
sign / (1 / std::abs(qoverpbrem) + std::abs(calomeots[i].energyLoss()->deltaE()));
4100 meff->setdelta_p(1000 * (qoverpbrem - qoverp));
4101 const AmgVector(5) & newpar = layerpar->parameters();
4103 prevtrackpars = layerpar->associatedSurface().createUniqueTrackParameters(
4104 newpar[0], newpar[1], newpar[2], newpar[3], qoverp, std::nullopt
4108 matstates.insert(matstates.begin(), std::make_unique<GXFTrackState>(std::move(meff), std::move(layerpar)));
4113 if (lasthit == lastmuonhit &&
cache.m_extmat) {
4114 std::unique_ptr<const Trk::TrackParameters> muonpar1;
4116 if (lastcalopar !=
nullptr) {
4118 if (msEntranceIsValid) {
4122 *
cache.m_msEntrance,
4126 if (muonpar1 !=
nullptr) {
4134 rot.col(2) = trackdir;
4136 trans.linear().matrix() << rot;
4137 trans.translation() << muonpar1->position() - .1 * trackdir;
4140 std::unique_ptr<const TrackParameters> curvlinpar(
m_extrapolator->extrapolateDirectly(
4148 if (curvlinpar !=
nullptr) {
4149 muonpar1 = std::move(curvlinpar);
4153 muonpar1 = std::unique_ptr<TrackParameters>(lastcalopar->
clone());
4157 muonpar1 = std::unique_ptr<TrackParameters>(refpar->clone());
4162 if (muonpar1 !=
nullptr) {
4164 muonpar1->position(),
4165 muonpar1->momentum().unit()
4169 double distance = getDistance(distsol);
4172 0) and (firstmuonhit !=
nullptr)) {
4174 muonpar1->position(),
4175 muonpar1->momentum().unit()
4181 distance = distsol.
first();
4184 std::abs(distsol.
first()) < std::abs(distsol.
second()) ?
4190 if (distance < 0 && distsol.
numberOfSolutions() > 0 && (firstidhit ==
nullptr)) {
4191 if (firstmuonpar !=
nullptr) {
4192 AmgVector(5) newpars = firstmuonpar->parameters();
4199 newpars[0], newpars[1], newpars[2], newpars[3], newpars[4], std::nullopt
4202 std::unique_ptr<const TrackParameters> tmppar(
m_propagator->propagateParameters(
4212 if (tmppar !=
nullptr) {
4213 muonpar1 = std::move(tmppar);
4219 ATH_MSG_DEBUG(
"Obtaining downstream layers from Extrapolator");
4220 if (matvec_used)
cache.m_matTempStore.push_back( std::move(matvec) );
4223 states.back()->associatedSurface(),
4227 matvec_used =
false;
4234 if (matvec && !matvec->empty()) {
4235 for (
int j = 0; j < (int) matvec->size(); j++) {
4241 std::unique_ptr<GXFMaterialEffects> meff = std::make_unique<GXFMaterialEffects>(*meot);
4246 std::abs(meff->deltaE()) > 25 &&
4247 std::abs((*matvec)[j]->trackParameters()->position().z()) < 13000
4252 const TrackParameters * newparams = (*matvec)[j]->trackParameters() !=
nullptr ? (*matvec)[j]->trackParameters()->clone() :
nullptr;
4254 matstates.push_back(std::make_unique<GXFTrackState>(
4256 std::unique_ptr<const TrackParameters>(newparams)
4266 if (firsthit == firstmuonhit &&
cache.m_extmat && (firstcalopar !=
nullptr)) {
4267 std::unique_ptr<const Trk::TrackParameters> muonpar1;
4270 if (msEntranceIsValid) {
4274 *
cache.m_msEntrance,
4278 if (muonpar1 !=
nullptr) {
4286 rot.col(2) = trackdir;
4288 trans.linear().matrix() << rot;
4289 trans.translation() << muonpar1->position() - .1 * trackdir;
4292 std::unique_ptr<const TrackParameters> curvlinpar(
m_extrapolator->extrapolateDirectly(
4300 if (curvlinpar !=
nullptr) {
4301 muonpar1 = std::move(curvlinpar);
4305 muonpar1 = std::unique_ptr<const TrackParameters>(firstcalopar->
clone());
4311 if (muonpar1 !=
nullptr) {
4313 muonpar1->position(),
4314 muonpar1->momentum().unit()
4318 const double distance = getDistance(distsol);
4322 ATH_MSG_DEBUG(
"Collecting upstream muon material from extrapolator");
4323 if (matvec_used)
cache.m_matTempStore.push_back( std::move(matvec) );
4326 states[0]->associatedSurface(),
4330 matvec_used =
false;
4332 if (matvec && !matvec->empty()) {
4333 ATH_MSG_DEBUG(
"Retrieved " << matvec->size() <<
" material states");
4335 for (
int j = 0; j < (int) matvec->size(); j++) {
4338 if (meb !=
nullptr) {
4343 std::unique_ptr<GXFMaterialEffects> meff = std::make_unique<GXFMaterialEffects>(*meot);
4348 std::abs(meff->deltaE()) > 25 &&
4349 std::abs((*matvec)[j]->trackParameters()->position().z()) < 13000
4355 (*matvec)[j]->trackParameters() !=
nullptr
4356 ? (*matvec)[j]->trackParameters()->clone()
4359 matstates.insert(matstates.begin(), std::make_unique<GXFTrackState>(
4361 std::unique_ptr<const TrackParameters>(tmpparams)
4374 std::vector<std::unique_ptr<GXFTrackState>> & newstates = states;
4375 std::vector<std::unique_ptr<GXFTrackState>> oldstates = std::move(newstates);
4378 newstates.reserve(oldstates.size() + matstates.size());
4381 int firstlayerno = -1;
4383 if (
cache.m_acceleration) {
4387 const double cosphi = std::cos(refpar->parameters()[
Trk::phi0]);
4388 const double sinphi = std::sin(refpar->parameters()[
Trk::phi0]);
4390 for (
int i =
cache.m_acceleration ? 1 : 0; i < (
int) oldstates.size(); i++) {
4391 bool addlayer =
true;
4393 while (addlayer && layerno < (
int) matstates.size()) {
4395 const TrackParameters *layerpar = matstates[layerno]->trackParameters();
4397 const DistanceSolution distsol = oldstates[i]->associatedSurface().straightLineDistanceEstimate(
4402 const double distance = getDistance(distsol);
4410 const double trackimpact = std::abs(-refpar->position().x() * sinphi + refpar->position().y() * cosphi);
4412 if (trackimpact > cylinderradius - 5 * mm) {
4418 if (i == (
int) oldstates.size() - 1) {
4427 ATH_MSG_DEBUG(
"Material state with excessive scattering, skipping it");
4438 if (firstlayerno < 0) {
4439 firstlayerno = layerno;
4444 if ((layerpar !=
nullptr) && matEffects !=
pion && matEffects !=
muon) {
4446 const Layer *lay =
nullptr;
4448 if (tvol !=
nullptr) {
4464 if (matvec_used)
cache.m_matTempStore.push_back( std::move(matvec) );
4477 if ((persurf !=
nullptr) && (!
cache.m_acceleration || persurf->
center().perp() > 5)) {
4478 const AmgVector(5) & pars = param.parameters();
4480 pars[0], pars[1], pars[2], pars[3], pars[4], std::nullopt
4484 if (
cache.m_acceleration) {
4489 std::unique_ptr<const TrackParameters> per(
m_extrapolator->extrapolate(
4490 Gaudi::Hive::currentContext(),param, tmppersf,
oppositeMomentum,
false, matEffects));
4492 if (per ==
nullptr) {
4493 ATH_MSG_DEBUG(
"Cannot make Perigee with starting parameters");
4497 if(std::abs(per->position().z())>5000.) {
4498 ATH_MSG_WARNING(
"Pathological perigee well outside of tracking detector!! Returning nullptr");
4506 const EventContext& ctx,
4513 ATH_MSG_DEBUG(
"--> entering GlobalChi2Fitter::myfit_helper");
4528 cache.m_lastiter = 0;
4541 if (trajectory.
nDOF() < 0) {
4546 cache.m_phiweight.clear();
4547 cache.m_firstmeasurement.clear();
4548 cache.m_lastmeasurement.clear();
4551 ATH_MSG_WARNING(
"Attempt to apply material corrections with q/p=0, reject track");
4565 std::unique_ptr<const TrackParameters> per =
makePerigee(
cache, param, matEffects);
4567 if (!
cache.m_acceleration && (per ==
nullptr)) {
4577 cache.m_acceleration &&
4581 ATH_MSG_WARNING(
"Tracking Geometry Service and/or Material Updator Tool not configured");
4584 cache.m_fastmat =
false;
4589 !
cache.m_acceleration ||
4592 addMaterial(ctx,
cache, trajectory, per !=
nullptr ? per.get() : ¶m, matEffects);
4595 ctx,
cache, trajectory, per !=
nullptr ? per.get() : ¶m, matEffects);
4607 for (
const auto & state : trajectory.
trackStates()) {
4611 state->materialEffects()->deltaE() == 0
4613 scatstate2 = state.get();
4618 scatstate = state.get();
4625 if ((scatstate ==
nullptr) || (scatstate2 ==
nullptr)) {
4626 throw std::logic_error(
"Invalid scatterer");
4631 const int nstates = (int) trajectory.
trackStates().size();
4639 std::unique_ptr<const TrackParameters> nearestpar;
4640 double mindist = 99999;
4641 std::vector < GXFTrackState * >mymatvec;
4644 if ((*it).trackParameters() ==
nullptr) {
4648 const double distance = persurf
4650 (*it).trackParameters()->position(),
4651 (*it).trackParameters()->momentum().unit())
4654 const bool insideid = (
4655 (
cache.m_caloEntrance ==
nullptr) ||
4656 cache.m_caloEntrance->inside((*it).trackParameters()->position())
4660 (((*it).measurement() !=
nullptr) && insideid) || (
4661 ((*it).materialEffects() !=
nullptr) &&
4663 (*it).materialEffects()->deltaE() == 0 ||
4664 ((*it).materialEffects()->sigmaDeltaPhi() == 0 &&
4666 (*it).materialEffects()->deltaPhi() != 0
4670 const double dist = ((*it).trackParameters()->
position() -
vertex).perp();
4671 if (dist < mindist) {
4679 if (((*it).materialEffects() !=
nullptr) && distance > 0) {
4680 mymatvec.push_back(it.get());
4684 if (nearestpar ==
nullptr) {
4688 for (
auto & state : mymatvec) {
4690 const Surface &matsurf = state->associatedSurface();
4692 nearestpar->position(), nearestpar->momentum().unit());
4694 const double distance = getDistance(distsol);
4700 std::unique_ptr<const TrackParameters> tmppar(
m_propagator->propagateParameters(
4710 if (tmppar ==
nullptr) {
4722 if (tmppar ==
nullptr) {
4732 AmgVector(5) newpars = tmppar->parameters();
4734 if (state->materialEffects()->sigmaDeltaE() > 0) {
4735 newpars[
Trk::qOverP] += .001 * state->materialEffects()->delta_p();
4738 const double de = std::abs(state->materialEffects()->deltaE());
4739 const double oldp = std::abs(1 / newpars[
Trk::qOverP]);
4740 const double newp2 = oldp * oldp - 2 * de * std::sqrt(mass * mass + oldp * oldp) + de * de;
4746 nearestpar = tmppar->associatedSurface().createUniqueTrackParameters(
4747 newpars[0], newpars[1], newpars[2], newpars[3], newpars[4], std::nullopt
4751 std::unique_ptr<Trk::TrackParameters> tmpPars(
m_propagator->propagateParameters(
4762 if (tmpPars !=
nullptr) {
4763 per.reset(
static_cast < const Perigee *
>(tmpPars.release()));
4768 const double oldp = 1. / std::abs(per->parameters()[
Trk::qOverP]);
4770 const double newp = std::sqrt(oldp * oldp + 2 * toteloss * std::sqrt(oldp * oldp + mass * mass) + toteloss * toteloss);
4771 AmgVector(5) params = per->parameters();
4774 per = per->associatedSurface().createUniqueTrackParameters(
4775 params[0], params[1], params[2], params[3], params[4], std::nullopt
4779 if (per ==
nullptr) {
4796 }
else if (per ==
nullptr) {
4813 const AmgVector(5) & pars = per->parameters();
4814 per = per->associatedSurface().createUniqueTrackParameters(
4815 pars[0], pars[1], pars[2], pars[3], 0, std::nullopt
4821 if (per !=
nullptr) {
4831 Eigen::MatrixXd a_inv;
4832 a.resize(nfitpar, nfitpar);
4837 derivPool.setZero();
4839 for (std::unique_ptr<GXFTrackState> & state : trajectory.
trackStates()) {
4840 if (state->materialEffects() !=
nullptr) {
4843 state->setDerivatives(derivPool);
4846 bool doderiv =
true;
4847 const int tmpminiter =
cache.m_miniter;
4849 for (
int it = 0; it <
m_maxit; ++it) {
4850 cache.m_lastiter = it;
4856 cache.m_miniter = tmpminiter;
4861 cache.m_fittercode =
4871 cache.m_miniter = tmpminiter;
4878 const double redchi2 = (trajectory.
nDOF() > 0) ? trajectory.
chi2() / trajectory.
nDOF() : 0;
4879 const double prevredchi2 = (trajectory.
nDOF() > 0) ? trajectory.
prevchi2() / trajectory.
nDOF() : 0;
4888 (redchi2 < prevredchi2 &&
4889 (redchi2 > prevredchi2 - 1 || redchi2 < 2)) ||
4890 nsihits + ntrthits == nhits
4895 if (it != 1 || nsihits != 0 || trajectory.
nDOF() <= 0 || trajectory.
chi2() / trajectory.
nDOF() <= 3) {
4900 cache.m_miniter = tmpminiter;
4910 if (ntrtprechits+ntrttubehits) {
4911 phf = float(ntrtprechits)/float(ntrtprechits+ntrttubehits);
4913 if (phf<m_minphfcut && it>=3) {
4914 if ((ntrtprechits+ntrttubehits)>=15) {
4918 ATH_MSG_DEBUG(
"Iter = " << it <<
" | nTRTStates = " << ntrthits
4919 <<
" | nTRTPrecHits = " << ntrtprechits
4920 <<
" | nTRTTubeHits = " << ntrttubehits
4921 <<
" | nOutliers = "
4930 cache.m_miniter = tmpminiter;
4939 cache.m_miniter = tmpminiter;
4941 if (trajectory.
prefit() == 0) {
4945 Eigen::LLT < Eigen::MatrixXd >
const lltOfW(
a);
4946 if (lltOfW.info() == Eigen::Success) {
4950 const int ncols =
a.cols();
4951 Amg::MatrixX const weightInvAMG = Amg::MatrixX::Identity(ncols, ncols);
4952 a_inv = lltOfW.solve(weightInvAMG);
4963 (runOutlier ||
cache.m_sirecal) &&
4971 if (traj != &trajectory) {
4976 finaltrajectory = traj;
4983 if (!
cache.m_acceleration && (finaltrajectory->
prefit() == 0)) {
4984 if (nperpars == 5) {
4985 for (
int i = 0; i <
a.cols(); i++) {
4986 a_inv(4, i) *= .001;
4987 a_inv(i, 4) *= .001;
4991 int scatterPos = nperpars + 2 * nscat;
4992 for (
int bremno = 0; bremno < nbrem; bremno++, scatterPos++) {
4993 for (
int i = 0; i <
a.cols(); i++) {
4994 a_inv(scatterPos, i) *= .001;
4995 a_inv(i, scatterPos) *= .001;
5002 for (
int i = 0; i < nperparams; i++) {
5003 for (
int j = 0; j < nperparams; j++) {
5004 (errmat) (j, i) = a_inv(j, i);
5009 (errmat) (4, 4) = 1e-20;
5013 std::unique_ptr<const TrackParameters> measper(
5015 perpars[0], perpars[1], perpars[2], perpars[3], perpars[4], std::move(errmat)
5021 cache.m_fullcovmat = a_inv;
5025 std::unique_ptr<Track> track =
nullptr;
5027 if (finaltrajectory->
prefit() > 0) {
5028 if (finaltrajectory != &trajectory) {
5030 delete finaltrajectory;
5049 (track !=
nullptr) && (
5050 track->fitQuality()->numberDoF() != 0 &&
5051 track->fitQuality()->chiSquared() / track->fitQuality()->numberDoF() > cut
5054 track.reset(
nullptr);
5058 if (track ==
nullptr) {
5062 if (finaltrajectory != &trajectory) {
5063 delete finaltrajectory;
5066 return track.release();
5070 const EventContext& ctx,
5075 int & bremno_maxbrempull,
5080 std::vector<std::unique_ptr<GXFTrackState>> & states = trajectory.
trackStates();
5100 const int nmeas = (int)
res.size();
5110 const int nDOF = trajectory.
nDOF();
5111 const bool doNewPseudoMeasurements = (
5115 std::abs((trajectory.
prevchi2() - trajectory.
chi2()) / nDOF) < 15 &&
5117 (nperpars == 0 || nidhits > 0)
5128 double maxbrempull = -0.2;
5137 for (
int hitno = 0; hitno < (int) states.size(); hitno++) {
5138 std::unique_ptr<GXFTrackState> & state = states[hitno];
5152 doNewPseudoMeasurements &&
5154 !state->associatedSurface().isFree() &&
5155 !state->isRecalibrated()
5158 covMatrix(0, 0) = 100;
5160 std::unique_ptr<const PseudoMeasurementOnTrack> newpseudo = std::make_unique<const PseudoMeasurementOnTrack>(
5162 std::move(covMatrix),
5166 state->setMeasurement(std::move(newpseudo));
5167 measbase = state->measurement();
5174 double *errors = state->measurementErrors();
5176 for (
int i = 0; i < 5; i++) {
5196 res[measno] = residuals[i];
5202 res[measno] = -std::remainder(-
res[measno], 2 *
M_PI);
5211 double *errors = state->measurementErrors();
5212 for (
int i = 0; i < 5; i++) {
5213 if (errors[i] > 0) {
5214 error[measno] = errors[i];
5225 ((trajectory.
prefit() == 0) || state->materialEffects()->deltaE() == 0)
5227 const double deltaPhi = state->materialEffects()->deltaPhi();
5228 const double measDeltaPhi = state->materialEffects()->measuredDeltaPhi();
5229 const double sigma2deltaPhi = std::pow(state->materialEffects()->sigmaDeltaPhi(), 2);
5230 const double deltaTheta = state->materialEffects()->deltaTheta();
5231 const double sigma2deltaTheta = std::pow(state->materialEffects()->sigmaDeltaTheta(), 2);
5233 if (trajectory.
prefit() != 1) {
5234 b[nperpars + 2 * scatno] -= (
deltaPhi - measDeltaPhi) / sigma2deltaPhi;
5235 b[nperpars + 2 * scatno + 1] -= deltaTheta / sigma2deltaTheta;
5237 b[nperpars + scatno] -= deltaTheta / sigma2deltaTheta;
5242 deltaTheta * deltaTheta / sigma2deltaTheta
5251 if ((state->materialEffects() !=
nullptr) && state->materialEffects()->sigmaDeltaE() > 0) {
5252 double averagenergyloss = std::abs(state->materialEffects()->deltaE());
5253 const double qoverpbrem = limitInversePValue(1000 * states[hitno]->trackParameters()->parameters()[
Trk::qOverP]);
5254 const double qoverp = limitInversePValue(qoverpbrem - state->materialEffects()->delta_p());
5255 const double pbrem = 1. / std::abs(qoverpbrem);
5256 const double p = 1. / std::abs(qoverp);
5257 const double mass = .001 * trajectory.
mass();
5258 const double energy = std::sqrt(p * p + mass * mass);
5259 const double bremEnergy = std::sqrt(pbrem * pbrem + mass * mass);
5261 const double resMaterial = .001 * averagenergyloss - energy + bremEnergy;
5262 res[nmeas - nbrem + bremno] = resMaterial;
5264 const double sigde = state->materialEffects()->sigmaDeltaE();
5265 const double sigdepos = state->materialEffects()->sigmaDeltaEPos();
5266 const double sigdeneg = state->materialEffects()->sigmaDeltaENeg();
5268 double errorMaterial = .001 * state->materialEffects()->sigmaDeltaE();
5269 error[nmeas - nbrem + bremno] = errorMaterial;
5280 if (state->materialEffects()->isKink()) {
5281 maxbrempull = -999999999;
5282 state_maxbrempull =
nullptr;
5286 cache.m_asymeloss &&
5288 trajectory.
prefit() == 0 &&
5290 sigde != sigdepos &&
5293 const double elosspull = resMaterial / errorMaterial;
5295 if (trajectory.
mass() > 100) {
5300 if (std::abs(elosspull) > 1) {
5301 if (elosspull < -1) {
5302 state->materialEffects()->setSigmaDeltaE(sigdepos);
5304 state->materialEffects()->setSigmaDeltaE(sigdeneg);
5307 errorMaterial = .001 * state->materialEffects()->sigmaDeltaE();
5308 error[nmeas - nbrem + bremno] = errorMaterial;
5320 !state->materialEffects()->isKink() && (
5321 (
m_fixbrem == -1 && elosspull < maxbrempull) ||
5325 bremno_maxbrempull = bremno;
5326 state_maxbrempull = state.get();
5327 maxbrempull = elosspull;
5336 trajectory.
prefit() == 0 &&
5337 state->materialEffects()->sigmaDeltaPhi() == 0 &&
5338 state->materialEffects()->isMeasuredEloss() &&
5339 resMaterial / (.001 * state->materialEffects()->sigmaDeltaEAve()) > 2.5
5341 const TrackParameters* parforcalo = states[hitno - 2]->trackParameters();
5344 std::vector<MaterialEffectsOnTrack> calomeots =
5357 if (calomeots.size() == 3) {
5358 averagenergyloss = std::abs(calomeots[1].energyLoss()->deltaE());
5359 const double newres = .001 * averagenergyloss - energy + bremEnergy;
5360 const double newerr = .001 * calomeots[1].energyLoss()->sigmaDeltaE();
5362 const double oldPull = resMaterial / errorMaterial;
5363 const double newPull = newres / newerr;
5365 if (std::abs(newPull) < std::abs(oldPull)) {
5366 ATH_MSG_DEBUG(
"Changing from measured to parametrized energy loss");
5368 state->materialEffects()->setEloss(std::unique_ptr<EnergyLoss>(calomeots[1].energyLoss()->clone()));
5369 state->materialEffects()->setSigmaDeltaE(calomeots[1].energyLoss()->sigmaDeltaE());
5370 res[nmeas - nbrem + bremno] = newres;
5371 error[nmeas - nbrem + bremno] = newerr;
5375 state->materialEffects()->setMeasuredEloss(
false);
5385 for (
int imeas = 0; imeas < nmeas; imeas++) {
5386 if (
error[imeas] == 0) {
5407 const double oldChi2 = trajectory.
prevchi2();
5408 const double newChi2 = trajectory.
chi2();
5413 const double nDOF = trajectory.
nDOF();
5414 const double oldRedChi2 = (nDOF > 0) ? oldChi2 / nDOF : 0;
5415 const double newRedChi2 = (nDOF > 0) ? newChi2 / nDOF : 0;
5418 trajectory.
prefit() > 0 && (
5419 (newRedChi2 < 2 && it != 0) ||
5420 (newRedChi2 < oldRedChi2 + .1 && std::abs(newRedChi2 - oldRedChi2) < 1 && it != 1)
5433 int miniter = (nsihits != 0 && nsihits + ntrthits == nhits) ? 1 : 2;
5434 miniter = std::max(miniter,
cache.m_miniter);
5436 if (it >= miniter && std::abs(oldChi2 - newChi2) < 1) {
5443 const int bremno_maxbrempull,
5449 if (state_maxbrempull ==
nullptr) {
5461 const int nmeas = (int)
res.size();
5464 const double oldError =
error[nmeas - nbrem + bremno_maxbrempull];
5466 error[nmeas - nbrem + bremno_maxbrempull] = newError;
5469 if (
a.cols() != nFitPars) {
5473 const double errorRatio = oldError / newError;
5474 const double errorReductionRatio = 1 - std::pow(errorRatio, 2);
5477 for (
int i = 0; i < nFitPars; i++) {
5478 if (weightderiv(nmeas - nbrem + bremno_maxbrempull, i) == 0) {
5482 for (
int j = i; j < nFitPars; j++) {
5483 const double newaij =
a(i, j) - errorReductionRatio *
5484 weightderiv(nmeas - nbrem + bremno_maxbrempull, i) *
5485 weightderiv(nmeas - nbrem + bremno_maxbrempull, j);
5487 a.fillSymmetric(i, j, newaij);
5489 weightderiv(nmeas - nbrem + bremno_maxbrempull, i) *= errorRatio;
5498 std::vector<std::unique_ptr<GXFTrackState>> & states = trajectory.
trackStates();
5511 const int nmeas = (int) weightderiv.rows();
5513 for (std::unique_ptr<GXFTrackState> & state : states) {
5517 const auto [scatmin, scatmax] = std::minmax(scatno, nscatupstream);
5518 const auto [bremmin, bremmax] = std::minmax(bremno, nbremupstream);
5525 const double sinStereo =
5527 state->sinStereo() :
5529 const double cosStereo =
5531 std::sqrt(1 - std::pow(sinStereo, 2)) :
5539 auto getThisDeriv = [sinStereo, cosStereo, &derivatives](
int i,
int j) ->
double {
5540 if (i == 0 && sinStereo != 0) {
5541 return derivatives(0, j) * cosStereo + sinStereo * derivatives(1, j);
5543 return derivatives(i, j);
5547 for (
int i = 0; i < 5; i++) {
5562 if (i == 0 && sinStereo != 0) {
5563 weightderiv.row(measno).head(cols) =
5564 (derivatives.row(0).
head(cols) * cosStereo +
5565 sinStereo * derivatives.row(1).head(cols)) /
5568 weightderiv.row(measno).head(cols) = derivatives.row(i).head(cols) /
error[measno];
5572 for (
int j = scatmin; j < scatmax; j++) {
5573 if (trajectory.
prefit() == 1) {
5574 const int index = nperparams + j;
5577 const int index = nperparams + 2 * j;
5579 weightderiv(measno,
index + 1) = getThisDeriv(i,
index + 1) /
error[measno];
5583 for (
int j = bremmin; j < bremmax; j++) {
5584 const int index = j + nperparams + 2 * nscat;
5591 double *errors = state->measurementErrors();
5592 for (
int i = 0; i < 5; i++) {
5593 if (errors[i] > 0) {
5599 ((trajectory.
prefit() == 0) || state->materialEffects()->deltaE() == 0)
5604 if ((state->materialEffects() !=
nullptr) && state->materialEffects()->sigmaDeltaE() > 0) {
5606 const double qoverpbrem = limitInversePValue(1000 * state->trackParameters()->parameters()[
Trk::qOverP]);
5607 const double qoverp = limitInversePValue(qoverpbrem - state->materialEffects()->delta_p());
5609 const double mass = .001 * trajectory.
mass();
5611 const auto thisMeasurementIdx{nmeas - nbrem + bremno};
5615 auto multiplier = [] (
double mass,
double qOverP){
5618 const auto qoverpTerm {multiplier(mass, qoverp) /
error[thisMeasurementIdx]};
5619 const auto qoverpBremTerm {multiplier(mass, qoverpbrem) /
error[thisMeasurementIdx]};
5622 weightderiv(thisMeasurementIdx, 4) = qoverpBremTerm - qoverpTerm;
5625 const auto bremNoBase = nperparams + 2 * nscat;
5626 if (bremno < nbremupstream) {
5627 weightderiv(thisMeasurementIdx, bremNoBase + bremno) = qoverpTerm;
5628 for (
int bremno2 = bremno + 1; bremno2 < nbremupstream; bremno2++) {
5629 weightderiv(thisMeasurementIdx, bremNoBase + bremno2) = qoverpTerm - qoverpBremTerm;
5632 weightderiv(thisMeasurementIdx, bremNoBase + bremno) = qoverpBremTerm;
5633 for (
int bremno2 = nbremupstream; bremno2 < bremno; bremno2++) {
5634 weightderiv(thisMeasurementIdx, bremNoBase + bremno2) = qoverpBremTerm - qoverpTerm;
5653 const int nMeas = (int)
res.size();
5655 cache.m_firstmeasurement.resize(nFitPars);
5656 cache.m_lastmeasurement.resize(nFitPars);
5658 for (
int i = 0; i < nPerPars; i++) {
5659 cache.m_firstmeasurement[i] = 0;
5660 cache.m_lastmeasurement[i] = nMeas - nBrem;
5667 const std::unique_ptr<GXFTrackState> & state = trajectory.
trackStates()[i];
5670 if (meff ==
nullptr) {
5671 measno += state->numberOfMeasuredParameters();
5675 const int firstMeasurement = i < nUpstreamStates ? 0 : measno;
5676 const int lastMeasurement = i < nUpstreamStates ? measno : nMeas - nBrem;
5679 && (trajectory.
prefit() == 0 || meff->
deltaE() == 0)) {
5680 const int scatterPos = nPerPars + 2 * scatno;
5682 cache.m_firstmeasurement[scatterPos] = firstMeasurement;
5683 cache.m_lastmeasurement[scatterPos] = lastMeasurement;
5685 cache.m_firstmeasurement[scatterPos + 1] = firstMeasurement;
5686 cache.m_lastmeasurement[scatterPos + 1] = lastMeasurement;
5692 const int bremPos = nPerPars + nScatPars + bremno;
5694 cache.m_firstmeasurement[bremPos] = firstMeasurement;
5695 cache.m_lastmeasurement[bremPos] = lastMeasurement;
5716 const int nMeas = (int)
res.size();
5718 for (
int k = 0; k < nFitPars; k++) {
5719 const int minMeasK =
cache.m_firstmeasurement[k];
5720 const int maxMeasK =
cache.m_lastmeasurement[k];
5727 for (
int measno = minMeasK; measno < maxMeasK; measno++) {
5728 b[k] +=
res[measno] * (1. /
error[measno]) * weightDeriv(measno, k);
5738 if (k == 4 || k >= nPerPars + nScatPars) {
5739 for (
int measno = nMeas - nBrem; measno < nMeas; measno++) {
5740 b[k] +=
res[measno] * (1. /
error[measno]) * weightDeriv(measno, k);
5754 for (
int k = 0; k < nFitPars; k++) {
5755 for (
int l = k; l < nFitPars; l++) {
5756 const int minMeas = std::max(
cache.m_firstmeasurement[k],
cache.m_firstmeasurement[l]);
5757 const int maxMeas = std::min(
cache.m_lastmeasurement[k],
cache.m_lastmeasurement[l]);
5760 for (
int measno = minMeas; measno < maxMeas; measno++) {
5761 a_kl += weightDeriv(measno, k) * weightDeriv(measno, l);
5764 a.fillSymmetric(l, k, a_kl);
5782 const int nMeas = (int)
res.size();
5789 for (
int k = nPerPars; k < nPerPars + nScatPars; k += 2) {
5790 a(k, k) += 1. / std::pow(scatSigmas[scatno].first, 2);
5791 a(k + 1, k + 1) += 1. / std::pow(scatSigmas[scatno].second, 2);
5799 for (
int measno = nMeas - nBrem; measno < nMeas; measno++) {
5800 for (
int k = 4; k < nFitPars; k++) {
5802 k = nPerPars + nScatPars;
5805 for (
int l = k; l < nFitPars; l++) {
5807 l = nPerPars + nScatPars;
5810 const double a_kl =
a(l, k) + weightDeriv(measno, k) * weightDeriv(measno, l);
5811 a.fillSymmetric(l, k, a_kl);
5823 const double oldRedChi2,
5824 const double newRedChi2
5832 bool weightChanged =
false;
5839 double newPhiWeight = 1.1;
5840 double newThetaWeight = 1.001;
5841 if (trajectory.
prefit() == 0) {
5847 newPhiWeight = 1.00000001;
5848 }
else if (it == 1) {
5849 newPhiWeight = 1.0000001;
5850 }
else if (it <= 3) {
5851 newPhiWeight = 1.0001;
5852 }
else if (it <= 6) {
5853 newPhiWeight = 1.01;
5856 if (newRedChi2 > oldRedChi2 - 1 && newRedChi2 < oldRedChi2) {
5857 newPhiWeight = 1.0001;
5858 newThetaWeight = 1.0001;
5859 }
else if (newRedChi2 > oldRedChi2 - 25 && newRedChi2 < oldRedChi2) {
5860 newPhiWeight = 1.001;
5861 newThetaWeight = 1.0001;
5868 std::size_t scatno = 0;
5873 for (
const auto & state : trajectory.
trackStates()) {
5876 if (meff ==
nullptr) {
5880 const bool isValidPlaneSurface =
5882 static_cast<const PlaneSurface *
>(&state->associatedSurface()) !=
nullptr;
5887 if (meff->
deltaE() == 0 || (trajectory.
prefit() == 0 && isValidPlaneSurface)) {
5888 weightChanged =
true;
5890 const int scatNoIndex = 2 * scatno + nPerPars;
5893 if (scatno >=
cache.m_phiweight.size()) {
5894 std::stringstream message;
5895 message <<
"scatno is out of range " << scatno <<
" !< " <<
cache.m_phiweight.size();
5896 throw std::range_error(message.str());
5904 a(scatNoIndex, scatNoIndex) /=
cache.m_phiweight[scatno];
5907 cache.m_phiweight[scatno] = newPhiWeight;
5908 a(scatNoIndex, scatNoIndex) *= newPhiWeight;
5909 }
else if (trajectory.
prefit() >= 2) {
5910 a(scatNoIndex, scatNoIndex) *= newPhiWeight;
5911 a(scatNoIndex + 1, scatNoIndex + 1) *= newThetaWeight;
5938 trajectory.
prefit() == 2 &&
5941 (newRedChi2 < oldRedChi2 - 25 || newRedChi2 > oldRedChi2)
5946 return weightChanged;
5955 std::size_t scatno = 0;
5964 if (scatno >=
cache.m_phiweight.size()) {
5965 std::stringstream message;
5966 message <<
"scatno is out of range " << scatno <<
" !< " <<
cache.m_phiweight.size();
5967 throw std::range_error(message.str());
5970 const bool isValidPlaneSurface =
5972 static_cast<const PlaneSurface *
>(&state->associatedSurface()) !=
nullptr;
5974 if (meff->
deltaE() == 0 || isValidPlaneSurface) {
5975 const int scatNoIndex = 2 * scatno + nPerPars;
5976 a(scatNoIndex, scatNoIndex) /=
cache.m_phiweight[scatno];
5977 cache.m_phiweight[scatno] = 1;
5991 const EventContext& ctx,
6000 const int nDOFold = trajectory.
nDOF();
6001 const double oldChi2 = trajectory.
chi2();
6002 const double oldRedChi2 = nDOFold > 0 ? oldChi2 / nDOFold : 0;
6004 if (
cache.m_phiweight.empty()) {
6024 int bremno_maxbrempull = 0;
6040 if ((state_maxbrempull !=
nullptr) && trajectory.
converged()) {
6049 const int nDOFnew = trajectory.
nDOF();
6050 const double newChi2 = trajectory.
chi2();
6051 const double newRedChi2 = nDOFnew > 0 ? newChi2 / nDOFnew : 0;
6053 ATH_MSG_DEBUG(
"old chi2: " << oldChi2 <<
"/" << nDOFold <<
"=" << oldRedChi2 <<
6054 ", new chi2: " << newChi2 <<
"/" << nDOFnew <<
"=" << newRedChi2);
6065 if (
cache.m_firstmeasurement.empty()) {
6089 if (doDeriv || weightChanged) {
6100 if (trajectory.
prefit() == 0) {
6106 if (nSiHits + nTrtHits !=
nHits) {
6113 (newRedChi2 < 2 || (newRedChi2 < oldRedChi2 && newRedChi2 > oldRedChi2 - .5))
6135 Eigen::LLT<Eigen::MatrixXd>
const llt(lu_m);
6137 if (llt.info() != Eigen::Success) {
6161 double d0 = refpar->parameters()[
Trk::d0];
6162 double z0 = refpar->parameters()[
Trk::z0];
6165 double qoverp = refpar->parameters()[
Trk::qOverP];
6167 if (nperparams > 0) {
6168 d0 += deltaParameters[0];
6169 z0 += deltaParameters[1];
6170 phi += deltaParameters[2];
6171 theta += deltaParameters[3];
6172 qoverp = (trajectory.
m_straightline) ? 0 : .001 * deltaParameters[4] + qoverp;
6184 std::vector < std::pair < double, double >>&scatangles = trajectory.
scatteringAngles();
6185 for (
int i = 0; i < nscat; i++) {
6186 scatangles[i].first += deltaParameters[2 * i + nperparams];
6187 scatangles[i].second += deltaParameters[2 * i + nperparams + 1];
6193 std::vector < double >&delta_ps = trajectory.
brems();
6194 for (
int i = 0; i < nbrem; i++) {
6195 delta_ps[i] += deltaParameters[nperparams + 2 * nscat + i];
6201 std::unique_ptr<const TrackParameters> newper(
6221 const EventContext& evtctx
6232 if (!splitProbContainer.
isValid()) {
6236 std::vector<std::unique_ptr<GXFTrackState>> & states = trajectory.
trackStates();
6243 for (
size_t stateno = 0; stateno < states.size(); stateno++) {
6248 measno += states[stateno-1]->numberOfMeasuredParameters();
6251 std::unique_ptr<GXFTrackState> & state = states[stateno];
6263 const auto *
const rot =
static_cast<const RIO_OnTrack *
>(pMeas);
6264 prd = rot->prepRawData();
6274 const auto &splitProb = splitProbContainer->splitProbability(
pixelCluster);
6275 if (!splitProb.isSplit()) {
6276 ATH_MSG_DEBUG(
"Pixel cluster is not split so no need to update" );
6280 std::unique_ptr < const RIO_OnTrack > newrot;
6281 double *olderror = state->measurementErrors();
6284 double newerror[5] = {-1,-1,-1,-1,-1};
6285 double newres[2] = {-1,-1};
6287 newrot.reset(
m_ROTcreator->correct(*prd, *trackpars, evtctx));
6292 const Amg::MatrixX & covmat = newrot->localCovariance();
6294 newerror[0] = std::sqrt(covmat(0, 0));
6295 newres[0] = newrot->localParameters()[
Trk::locX] - trackpars->parameters()[
Trk::locX];
6296 newerror[1] = std::sqrt(covmat(1, 1));
6297 newres[1] = newrot->localParameters()[
Trk::locY] - trackpars->parameters()[
Trk::locY];
6299 if (
a.cols() != nfitpars) {
6304 for(
int k =0; k<2; k++ ){
6305 const double oldres =
res[measno+k];
6306 res[measno+k] = newres[k];
6307 err[measno+k] = newerror[k];
6309 for (
int i = 0; i < nfitpars; i++) {
6310 if (weightderiv(measno+k, i) == 0) {
6314 b[i] -= weightderiv(measno+k, i) * (oldres / olderror[k] - (newres[k] * olderror[k]) / (newerror[k] * newerror[k]));
6316 for (
int j = i; j < nfitpars; j++) {
6320 weightderiv(measno+k, i) *
6321 weightderiv(measno+k, j) *
6322 ((olderror[k] * olderror[k]) / (newerror[k] * newerror[k]) - 1)
6326 weightderiv(measno+k, i) *= olderror[k] / newerror[k];
6330 state->setMeasurement(std::move(newrot));
6331 state->setMeasurementErrors(newerror);
6346 const EventContext& ctx
6354 std::vector<std::unique_ptr<GXFTrackState>> & states = trajectory.
trackStates();
6360 if (
a.cols() != nfitpars) {
6368 bool outlierremoved =
false;
6369 bool hitrecalibrated =
false;
6371 for (
int stateno = 0; stateno < (int) states.size(); stateno++) {
6372 std::unique_ptr<GXFTrackState> & state = states[stateno];
6380 std::abs(state->trackParameters()->parameters()[
Trk::driftRadius]) > 1.05 * state->associatedSurface().bounds().r()
6385 outlierremoved =
true;
6387 double *errors = state->measurementErrors();
6388 const double olderror = errors[0];
6392 for (
int i = 0; i < nfitpars; i++) {
6393 if (weightderiv(measno, i) == 0) {
6397 b[i] -=
res[measno] * weightderiv(measno, i) / olderror;
6399 for (
int j = i; j < nfitpars; j++) {
6402 a(i, j) - weightderiv(measno, i) * weightderiv(measno, j)
6405 weightderiv(measno, i) = 0;
6409 }
else if (trtrecal) {
6410 double *errors = state->measurementErrors();
6411 const double olderror = errors[0];
6413 const auto *
const thisMeasurement{state->measurement()};
6422 const double trackradius = state->trackParameters()->parameters()[
Trk::driftRadius];
6424 std::unique_ptr<const Trk::RIO_OnTrack> newrot =
nullptr;
6425 const double distance = std::abs(std::abs(trackradius) - dcradius);
6427 if (distance < scalefactor * dcerror && (olderror > 1. || trackradius * oldradius < 0)) {
6429 }
else if (distance > scalefactor * dcerror && olderror < 1.) {
6433 if (newrot !=
nullptr) {
6435 hitrecalibrated =
true;
6439 if ((measno < 0) or (measno >= (
int)
res.size())) {
6440 throw std::runtime_error(
6441 "'res' array index out of range in TrkGlobalChi2Fitter/src/GlobalChi2Fitter.cxx:" + std::to_string(__LINE__)
6445 const double oldres =
res[measno];
6446 const double newres = newradius - state->trackParameters()->parameters()[
Trk::driftRadius];
6447 errors[0] = newerror;
6448 state->setMeasurement(std::move(newrot));
6452 for (
int i = 0; i < nfitpars; i++) {
6453 if (weightderiv(measno, i) == 0) {
6457 b[i] -= weightderiv(measno, i) * (oldres / olderror - (newres * olderror) / (newerror * newerror));
6459 for (
int j = i; j < nfitpars; j++) {
6463 !
cache.m_phiweight.empty() &&
6466 i < nperpars + 2 * nscats &&
6467 (i - nperpars) % 2 == 0
6469 weight =
cache.m_phiweight[(i - nperpars) / 2];
6474 a(i, j) + weightderiv(measno, i) * weightderiv(measno, j) * ((olderror * olderror) / (newerror * newerror) - 1) * weight
6477 weightderiv(measno, i) *= olderror / newerror;
6480 res[measno] = newres;
6481 err[measno] = newerror;
6490 measno += state->numberOfMeasuredParameters();
6494 if (trajectory.
nDOF() < 0) {
6499 if (outlierremoved || hitrecalibrated) {
6503 cache.m_miniter = it + 2;
6508 const EventContext& ctx,
6516 bool trackok =
false;
6518 std::unique_ptr < GXFTrajectory > cleanup_oldtrajectory;
6520 std::unique_ptr < GXFTrajectory > cleanup_newtrajectory;
6526 while (!trackok && oldtrajectory->
nDOF() > 0) {
6528 std::vector<std::unique_ptr<GXFTrackState>> & states = oldtrajectory->
trackStates();
6536 if (nhits != nsihits) {
6540 double maxsipull = -1;
6542 int hitno_maxsipull = -1;
6543 int measno_maxsipull = -1;
6544 int stateno_maxsipull = 0;
6556 for (
int stateno = 0; stateno < (int) states.size(); stateno++) {
6557 std::unique_ptr<GXFTrackState> & state = states[stateno];
6563 double *errors = state->measurementErrors();
6565 const Amg::MatrixX & hitcov = state->measurement()->localCovariance();
6566 const double sinstereo = state->sinStereo();
6567 const double cosstereo = (sinstereo == 0) ? 1 : std::sqrt(1 - sinstereo * sinstereo);
6568 double weight1 = -1;
6570 if (hitcov(0, 0) > trackcov(0, 0)) {
6571 if (sinstereo == 0) {
6572 weight1 = errors[0] * errors[0] - trackcov(0, 0);
6574 weight1 = errors[0] * errors[0] - (
6575 trackcov(0, 0) * cosstereo * cosstereo + 2 *
6576 trackcov(1, 0) * cosstereo * sinstereo + trackcov(1, 1) * sinstereo * sinstereo
6581 const double weight2 = (
6583 errors[1] * errors[1] - trackcov(1, 1) :
6587 double sipull1 = weight1 > 0 ? std::abs(
res[measno] / std::sqrt(weight1)) : -1;
6588 const double sipull2 = (
6590 std::abs(
res[measno + 1] / std::sqrt(weight2)) :
6593 sipull1 = std::max(sipull1, sipull2);
6595 if (sipull1 > maxsipull) {
6596 maxsipull = sipull1;
6597 measno_maxsipull = measno;
6598 state_maxsipull = state.get();
6599 stateno_maxsipull = stateno;
6600 hitno_maxsipull = hitno;
6611 measno += state->numberOfMeasuredParameters();
6615 const double maxpull = maxsipull;
6617 ATH_MSG_DEBUG(
" maxsipull: " << maxsipull <<
" hitno_maxsipull: " <<
6618 hitno_maxsipull <<
" n3sigma: " << n3sigma <<
" cut: " << cut <<
" cut2: " << cut2);
6629 state_maxsipull = oldtrajectory->
trackStates()[stateno_maxsipull].get();
6632 const auto *
const rot =
static_cast<const RIO_OnTrack *
>(pMeas);
6633 prd = rot->prepRawData();
6635 std::unique_ptr < const RIO_OnTrack > broadrot;
6640 Amg::VectorX parameterVector = trackpar_maxsipull->parameters();
6641 const std::unique_ptr<const TrackParameters> trackparForCorrect(
6650 state_maxsipull->trackCovariance())
6654 newerror[0] = newerror[1] = newerror[2] = newerror[3] = newerror[4] = -1;
6655 double newpull = -1;
6656 double newpull1 = -1;
6657 double newpull2 = -1;
6658 double newres1 = -1;
6659 double newres2 = -1;
6660 double newsinstereo = 0;
6673 const Amg::MatrixX & covmat = broadrot->localCovariance();
6675 if (state_maxsipull->
sinStereo() != 0) {
6676 const auto [covEigenValueSmall, covStereoAngle] = principalComponentAnalysis2x2(covmat);
6677 newerror[0] = std::sqrt(covEigenValueSmall);
6678 newsinstereo = std::sin(covStereoAngle);
6680 newerror[0] = std::sqrt(covmat(0, 0));
6683 const double cosstereo = (newsinstereo == 0) ? 1. : std::sqrt(1 - newsinstereo * newsinstereo);
6685 if (cosstereo != 1.) {
6687 cosstereo * (broadrot->localParameters()[
Trk::locX] - trackpar_maxsipull->parameters()[
Trk::locX]) +
6688 newsinstereo * (broadrot->localParameters()[
Trk::locY] - trackpar_maxsipull->parameters()[
Trk::locY])
6691 newres1 = broadrot->localParameters()[
Trk::locX] - trackpar_maxsipull->parameters()[
Trk::locX];
6694 if (newerror[0] == 0.0) {
6695 ATH_MSG_WARNING(
"Measurement error is zero or negative, treating as outlier");
6698 newpull1 = std::abs(newres1 / newerror[0]);
6702 newerror[1] = std::sqrt(covmat(1, 1));
6703 newres2 = broadrot->localParameters()[
Trk::locY] - trackpar_maxsipull->parameters()[
Trk::locY];
6704 newpull2 = std::abs(newres2 / newerror[1]);
6707 newpull = std::max(newpull1, newpull2);
6713 (newerror[0] > 1.5 * olderror[0] || newerror[1] > 1.5 * std::abs(olderror[1]))
6715 if ((measno_maxsipull < 0) or(measno_maxsipull >= (
int)
res.size())) {
6716 throw std::runtime_error(
6717 "'res' array index out of range in TrkGlobalChi2Fitter/src/GlobalChi2Fitter.cxx:" + std::to_string(__LINE__)
6722 newtrajectory = oldtrajectory;
6724 if (
a.cols() != nfitpars) {
6728 const double oldres1 =
res[measno_maxsipull];
6729 res[measno_maxsipull] = newres1;
6730 err[measno_maxsipull] = newerror[0];
6732 for (
int i = 0; i < nfitpars; i++) {
6733 if (weightderiv(measno_maxsipull, i) == 0) {
6737 b[i] -= weightderiv(measno_maxsipull, i) * (oldres1 / olderror[0] - (newres1 * olderror[0]) / (newerror[0] * newerror[0]));
6739 for (
int j = i; j < nfitpars; j++) {
6743 weightderiv(measno_maxsipull, i) *
6744 weightderiv(measno_maxsipull, j) *
6745 ((olderror[0] * olderror[0]) / (newerror[0] * newerror[0]) - 1)
6749 weightderiv(measno_maxsipull, i) *= olderror[0] / newerror[0];
6753 const double oldres2 =
res[measno_maxsipull + 1];
6754 res[measno_maxsipull + 1] = newres2;
6755 err[measno_maxsipull + 1] = newerror[1];
6757 for (
int i = 0; i < nfitpars; i++) {
6758 if (weightderiv(measno_maxsipull + 1, i) == 0) {
6762 b[i] -= weightderiv(measno_maxsipull + 1, i) * (oldres2 / olderror[1] - (newres2 * olderror[1]) / (newerror[1] * newerror[1]));
6764 for (
int j = i; j < nfitpars; j++) {
6768 weightderiv(measno_maxsipull + 1, i) *
6769 weightderiv(measno_maxsipull + 1, j) *
6770 ((olderror[1] * olderror[1]) / (newerror[1] * newerror[1]) - 1)
6775 weightderiv(measno_maxsipull + 1, i) *= olderror[1] / newerror[1];
6780 "Recovering outlier, hitno=" << hitno_maxsipull <<
" measno=" <<
6781 measno_maxsipull <<
" pull=" << maxsipull <<
" olderror_0=" <<
6782 olderror[0] <<
" newerror_0=" << newerror[0] <<
" olderror_1=" <<
6783 olderror[1] <<
" newerror_1=" << newerror[1]
6792 ((n3sigma < 2 && maxsipull > cut2 && maxsipull < cut) || n3sigma > 1) &&
6802 "Removing outlier, hitno=" << hitno_maxsipull <<
", measno=" <<
6803 measno_maxsipull <<
" pull=" << maxsipull
6810 cleanup_newtrajectory = std::make_unique<GXFTrajectory>(*oldtrajectory);
6811 newtrajectory = cleanup_newtrajectory.get();
6813 if (newa.cols() != nfitpars) {
6819 if ((measno_maxsipull < 0) or(measno_maxsipull >= (
int)
res.size())) {
6820 throw std::runtime_error(
6821 "'res' array index out of range in TrkGlobalChi2Fitter/src/GlobalChi2Fitter.cxx:" + std::to_string(__LINE__)
6825 const double oldres1 =
res[measno_maxsipull];
6826 newres[measno_maxsipull] = 0;
6828 for (
int i = 0; i < nfitpars; i++) {
6829 if (weightderiv(measno_maxsipull, i) == 0) {
6833 newb[i] -= weightderiv(measno_maxsipull, i) * oldres1 / olderror[0];
6835 for (
int j = i; j < nfitpars; j++) {
6839 weightderiv(measno_maxsipull, i) *
6840 weightderiv(measno_maxsipull, j)
6844 newweightderiv(measno_maxsipull, i) = 0;
6848 const double oldres2 =
res[measno_maxsipull + 1];
6849 newres[measno_maxsipull + 1] = 0;
6851 for (
int i = 0; i < nfitpars; i++) {
6852 if (weightderiv(measno_maxsipull + 1, i) == 0) {
6856 newb[i] -= weightderiv(measno_maxsipull + 1, i) * oldres2 / olderror[1];
6858 for (
int j = i; j < nfitpars; j++) {
6859 if (weightderiv(measno_maxsipull + 1, j) == 0) {
6866 weightderiv(measno_maxsipull + 1, i) *
6867 weightderiv(measno_maxsipull + 1, j)
6871 newweightderiv(measno_maxsipull + 1, i) = 0;
6875 newtrajectory->
setOutlier(stateno_maxsipull);
6889 for (
int it = 0; it <
m_maxit; ++it) {
6899 ctx,
cache, *newtrajectory, it, *newap, *newbp, lu_m, doderiv);
6915 const double oldchi2 = oldtrajectory->
chi2() / oldtrajectory->
nDOF();
6916 const double newchi2 = (newtrajectory->
nDOF() > 0) ? newtrajectory->
chi2() / newtrajectory->
nDOF() : 0;
6919 if (newtrajectory->
nDOF() != oldtrajectory->
nDOF() && maxsipull > cut2) {
6920 mindiff = (oldchi2 > .33 *
m_chi2cut || noutl > 0) ? .8 : 1.;
6922 if (noutl == 0 && maxsipull < cut - .5 && oldchi2 < .5 *
m_chi2cut) {
6927 if (newchi2 > oldchi2 || (newchi2 > oldchi2 - mindiff && newchi2 > .33 * oldchi2)) {
6928 ATH_MSG_DEBUG(
"Outlier not confirmed, keeping old trajectory");
6936 (void)cleanup_oldtrajectory.release();
6937 return oldtrajectory;
6939 if (oldtrajectory != newtrajectory) {
6940 cleanup_oldtrajectory = std::move(cleanup_newtrajectory);
6941 oldtrajectory = newtrajectory;
6948 Eigen::LLT < Eigen::MatrixXd >
const lltOfW(
a);
6949 if (lltOfW.info() == Eigen::Success) {
6953 const int ncols =
a.cols();
6954 Amg::MatrixX const weightInvAMG = Amg::MatrixX::Identity(ncols, ncols);
6955 fullcov = lltOfW.solve(weightInvAMG);
6973 oldtrajectory->
nDOF() > 0 &&
6982 (void)cleanup_oldtrajectory.release();
6983 return oldtrajectory;
6995 if (
const auto *pMeas{hit->measurement()};
7001 nrealmeas += hit->numberOfMeasuredParameters();
7005 cache.m_derivmat.setZero();
7011 if (
const auto *pMeas{hit->measurement()};
7017 for (
int i = measindex; i < measindex + hit->numberOfMeasuredParameters(); i++) {
7019 cache.m_derivmat(i, j) = derivs(measindex2, j) * errors[measindex2];
7020 if ((j == 4 && !oldtrajectory.
m_straightline) || j >= nperpars + 2 * nscat) {
7021 cache.m_derivmat(i, j) *= 1000;
7028 measindex += hit->numberOfMeasuredParameters();
7029 }
else if (hit->materialEffects() ==
nullptr) {
7030 measindex2 += hit->numberOfMeasuredParameters();
7036 const EventContext & ctx,
7044 std::unique_ptr<const TrackParameters> per(
nullptr);
7047 std::unique_ptr<const TrackParameters> prevpar(
7052 std::vector<std::pair<const Layer *, const Layer *>> & upstreamlayers = oldtrajectory.
upstreamMaterialLayers();
7055 for (
const auto & [layer1, layer2] : upstreamlayers | std::views::reverse) {
7056 if (prevpar ==
nullptr) {
7061 const Layer *layer = layer1 !=
nullptr ? layer1 : layer2;
7063 const DistanceSolution distsol = layer->surfaceRepresentation().straightLineDistanceEstimate(
7064 prevpar->position(), prevpar->momentum().unit()
7066 const double distance = getDistance(distsol);
7069 if (std::abs(distance) < 0.01) {
7073 if (distsol.
first() * distsol.
second() < 0 && !first) {
7078 if (first && distance > 0) {
7082 std::unique_ptr<const TrackParameters> layerpar(
7086 layer->surfaceRepresentation(),
7094 if (layerpar ==
nullptr) {
7098 if (layer->surfaceRepresentation().bounds().inside(layerpar->localPosition())) {
7102 prevpar = std::move(layerpar);
7115 if (startfactor > 0.5) {
7116 std::unique_ptr<const TrackParameters> updatedpar =
m_matupdator->update(
7120 if (updatedpar !=
nullptr) {
7140 if (endfactor > 0.5) {
7141 std::unique_ptr<const TrackParameters> updatedpar =
m_matupdator->update(
7145 if (updatedpar !=
nullptr) {
7151 if (prevpar !=
nullptr) {
7163 if (per ==
nullptr) {
7164 ATH_MSG_DEBUG(
"Failed to extrapolate to perigee, returning 0");
7183 std::unique_ptr<GXFTrackState>
7185 const EventContext & ctx,
7192 if (per ==
nullptr) {
7196 ATH_MSG_DEBUG(
"Final perigee: " << *per <<
" pos: " << per->position() <<
" pT: " << per->pT());
7202 const std::vector<std::unique_ptr<TrackParameters>> & hc,
7203 std::set<Identifier> & id_set,
7204 std::set<Identifier> & sct_set,
7214 for (
const std::unique_ptr<TrackParameters> & tp : hc) {
7220 if (tp ==
nullptr) {
7229 const TrkDetElementBase * de = tp->associatedSurface().associatedDetectorElement();
7231 if (de ==
nullptr) {
7242 if (id_set.find(
id) != id_set.end()) {
7261 }
else if (
m_DetID->is_sct(
id)) {
7272 }
else if (
m_DetID->is_sct(
id)) {
7281 const Identifier os = e->otherSide()->identify();
7291 if (sct_set.find(os) != sct_set.end()) {
7292 ++rv.m_sct_double_hole;
7321 for (
const std::unique_ptr<GXFTrackState> & s : trajectory.
trackStates()) {
7333 std::vector<std::reference_wrapper<GXFTrackState>> rv;
7341 for (
const std::unique_ptr<GXFTrackState> & s : trajectory.
trackStates()) {
7355 rv.emplace_back(*s);
7362 const TrkDetElementBase * de = s->trackParameters()->associatedSurface().associatedDetectorElement();
7364 if (de !=
nullptr) {
7377 if (s.get() == lastmeas) {
7387 const EventContext & ctx,
7388 const std::vector<std::reference_wrapper<GXFTrackState>> & states
7400 constexpr uint min_meas = 3;
7401 if (std::count_if(states.begin(), states.end(), [](
const GXFTrackState & s){ return s.getStateType(TrackStateOnSurface::Measurement); }) < min_meas) {
7405 bool seen_meas =
false;
7407 std::set<Identifier> id_set;
7408 std::set<Identifier> sct_set;
7414 for (std::size_t i = 0; i < states.size() - 1; i++) {
7437 const double dist = (beg.trackParameters()->
position() - end.trackParameters()->position()).norm();
7439 const bool zStartValid = std::abs(beg.trackParameters()->position().z())<10000.;
7441 ATH_MSG_DEBUG(
"Pathological track parameter well outside of detector");
7442 ATH_MSG_DEBUG(
"Propagator might have issue with this, skipping");
7443 ATH_MSG_VERBOSE(
"dumping track parameters " << *(beg.trackParameters()));
7452 if (seen_meas && dist >= 2.5 && zStartValid) {
7459 std::optional<std::vector<std::unique_ptr<TrackParameters>>> & hc = beg.getHoles();
7460 std::vector<std::unique_ptr<TrackParameters>> states;
7468 if (hc.has_value()) {
7469 states = std::move(*hc);
7507 std::vector<std::unique_ptr<Trk::TrackParameters>>
const bl =
m_extrapolator->extrapolateBlindly(
7528 const EventContext & ctx,
7534 auto trajectory = std::make_unique<Trk::TrackStates>();
7544 if (perigee_ts ==
nullptr) {
7550 trajectory->reserve(tmptrajectory.
trackStates().size());
7556 hit->resetTrackCovariance();
7562 hit->materialEffects()))
7568 auto trackState = hit->trackStateOnSurface();
7569 hit->resetTrackCovariance();
7570 trajectory->emplace_back(trackState.release());
7573 auto qual = std::make_unique<FitQuality>(tmptrajectory.
chi2(), tmptrajectory.
nDOF());
7593 std::unique_ptr<Track> rv = std::make_unique<Track>(info, std::move(trajectory), std::move(qual));
7603 std::unique_ptr<TrackSummary>
ts = std::make_unique<TrackSummary>();
7612 std::optional<TrackHoleCount> hole_count;
7619 std::vector<std::reference_wrapper<GXFTrackState>>
const states =
holeSearchStates(tmptrajectory);
7637 if (hole_count.has_value()) {
7650 rv->setTrackSummary(std::move(
ts));
7659 const EventContext & ctx,
7668 std::vector<std::unique_ptr<TrackParameters>> rv =
m_extrapolator->extrapolateStepwise(
7682 &rv.front()->associatedSurface() == &src.associatedSurface() ||
7683 trackParametersClose(*rv.front(), src, 0.001) ||
7687 rv.front().reset(
nullptr);
7698 &rv.back()->associatedSurface() == &src.associatedSurface() ||
7699 trackParametersClose(*rv.back(), src, 0.001) ||
7703 rv.back().reset(
nullptr);
7710 const EventContext & ctx,
7718 std::unique_ptr<const TrackParameters> rv;
7719 std::optional<TransportJacobian> jac{};
7730 if (rv !=
nullptr && calcderiv) {
7735 std::optional<std::vector<std::unique_ptr<TrackParameters>>> extrapolation;
7744 std::move(extrapolation)
7749 const EventContext & ctx,
7760 ctx, prev,
ts, propdir, bf, calcderiv, holesearch
7763 if (rv.m_parameters ==
nullptr) {
7764 propdir = invertPropdir(propdir);
7767 ctx, prev,
ts, propdir, bf, calcderiv, holesearch
7775 const EventContext& ctx,
7782 std::vector<std::unique_ptr<GXFTrackState>> & states = trajectory.
trackStates();
7785 std::unique_ptr<const TrackParameters> tmptrackpar;
7787 for (
int hitno = nstatesupstream - 1; hitno >= 0; hitno--) {
7788 const Surface &surf1 = states[hitno]->associatedSurface();
7795 const double distance = getDistance(distsol);
7817 (rv.m_parameters !=
nullptr) &&
7818 (prevtrackpar->
position() - rv.m_parameters->position()).mag() > 5 * mm
7824 if (rv.m_parameters ==
nullptr) {
7825 ATH_MSG_DEBUG(
"propagation failed, prev par: " << *prevtrackpar <<
7826 " pos: " << prevtrackpar->
position() <<
" destination surface: " << surf1);
7830 states[hitno]->setTrackParameters(std::move(rv.m_parameters));
7831 const TrackParameters *currenttrackpar = states[hitno]->trackParameters();
7832 const Surface &surf = states[hitno]->associatedSurface();
7834 if (rv.m_jacobian != std::nullopt) {
7836 states[hitno]->materialEffects() !=
nullptr &&
7837 states[hitno]->materialEffects()->deltaE() != 0 &&
7838 states[hitno]->materialEffects()->sigmaDeltaE() <= 0 &&
7841 const double p = 1. / std::abs(currenttrackpar->parameters()[
Trk::qOverP]);
7842 const double de = std::abs(states[hitno]->materialEffects()->deltaE());
7843 const double mass = trajectory.
mass();
7844 const double newp = std::sqrt(p * p + 2 * de * std::sqrt(mass * mass + p * p) + de * de);
7845 (*rv.m_jacobian) (4, 4) = ((p + p * de / std::sqrt(p * p + mass * mass)) / newp) * p * p / (newp * newp);
7848 states[hitno]->setJacobian(*rv.m_jacobian);
7849 }
else if (calcderiv) {
7856 if (meff !=
nullptr && hitno != 0) {
7858 surf, *meff, *states[hitno]->trackParameters(), trajectory.
mass(), -1
7861 if (std::holds_alternative<FitterStatusCode>(
r)) {
7862 return std::get<FitterStatusCode>(
r);
7865 tmptrackpar = std::move(std::get<std::unique_ptr<const TrackParameters>>(
r));
7866 prevtrackpar = tmptrackpar.get();
7868 prevtrackpar = currenttrackpar;
7874 for (
int hitno = nstatesupstream; hitno < (int) states.size(); hitno++) {
7875 const Surface &surf = states[hitno]->associatedSurface();
7879 const double distance = getDistance(distsol);
7896 (rv.m_parameters !=
nullptr) &&
7898 (prevtrackpar->
position() - rv.m_parameters->position()).mag() > 5 * mm
7903 if (rv.m_parameters ==
nullptr) {
7904 ATH_MSG_DEBUG(
"propagation failed, prev par: " << *prevtrackpar <<
7905 " pos: " << prevtrackpar->
7906 position() <<
" destination surface: " << surf);
7910 if (rv.m_jacobian != std::nullopt) {
7912 states[hitno]->materialEffects() !=
nullptr &&
7913 states[hitno]->materialEffects()->deltaE() != 0 &&
7914 states[hitno]->materialEffects()->sigmaDeltaE() <= 0 &&
7917 const double p = 1 / std::abs(rv.m_parameters->parameters()[
Trk::qOverP]);
7918 const double de = std::abs(states[hitno]->materialEffects()->deltaE());
7919 const double mass = trajectory.
mass();
7920 double newp = p * p - 2 * de * std::sqrt(mass * mass + p * p) + de * de;
7923 newp = std::sqrt(newp);
7926 (*rv.m_jacobian) (4, 4) = ((p - p * de / std::sqrt(p * p + mass * mass)) / newp) * p * p / (newp * newp);
7929 states[hitno]->setJacobian(*rv.m_jacobian);
7930 }
else if (calcderiv) {
7937 if (meff !=
nullptr) {
7939 surf, *meff, *rv.m_parameters, trajectory.
mass(), +1
7942 if (std::holds_alternative<FitterStatusCode>(
r)) {
7943 return std::get<FitterStatusCode>(
r);
7946 rv.m_parameters = std::move(std::get<std::unique_ptr<const TrackParameters>>(
r));
7949 states[hitno]->setTrackParameters(std::move(rv.m_parameters));
7950 prevtrackpar = states[hitno]->trackParameters();
7963 const AmgVector(5) & old = param.parameters();
7969 ATH_MSG_DEBUG(
"Angles out of range, phi: " << newphi <<
" theta: " << newtheta);
7973 double newqoverp = 0;
7979 const double oldp = std::abs(1 / old[
Trk::qOverP]);
7980 const double newp2 = oldp * oldp -
sign * 2 * std::abs(meff.
deltaE()) * std::sqrt(mass * mass + oldp * oldp) + meff.
deltaE() * meff.
deltaE();
7987 newqoverp = std::copysign(1 / std::sqrt(newp2), old[
Trk::qOverP]);
7994 old[0], old[1], newphi, newtheta, newqoverp, std::nullopt
8006 using Matrix55 = Eigen::Matrix<double, 5, 5>;
8008 Matrix55 initialjac;
8009 initialjac.setZero();
8010 initialjac(4, 4) = 1;
8012 Matrix55 jacvertex(initialjac);
8014 std::vector<Matrix55, Eigen::aligned_allocator<Matrix55>> jacscat(trajectory.
numberOfScatterers(), initialjac);
8015 std::vector<Matrix55, Eigen::aligned_allocator<Matrix55>> jacbrem(trajectory.
numberOfBrems(), initialjac);
8017 std::vector<std::unique_ptr<GXFTrackState>> & states = trajectory.
trackStates();
8026 for (
const bool forward : {
false,
true}) {
8028 hit_begin = nstatesupstream;
8029 hit_end = (int) states.size();
8030 scatno = nscatupstream;
8031 bremno = nbremupstream;
8033 hit_begin = nstatesupstream - 1;
8040 int hitno = hit_begin;
8041 forward ? (hitno < hit_end) : (hitno >= hit_end);
8042 hitno += (forward ? 1 : -1)
8045 state = states[hitno].get();
8049 if (fillderivmat && state->
derivatives().cols() != nfitpars) {
8057 const int jmaxbrem = 4;
8059 if (hitno == (forward ? hit_end - 1 : 0)) {
8060 if (!fillderivmat) {
8068 Eigen::Matrix<double, 5, 5> & jac = state->
jacobian();
8070 if (hitno == nstatesupstream + (forward ? 0 : -1)) {
8071 jacvertex.block<4, 5>(0, 0) = jac.block<4, 5>(0, 0);
8072 jacvertex(4, 4) = jac(4, 4);
8082 jcnt = jmax - jmin + 1;
8084 lp_bgn = forward ? nscatupstream : nscatupstream - 1;
8087 for (
int i = lp_bgn; forward ? (i < lp_end) : (i > lp_end); i += (forward ? 1 : -1)) {
8089 i == scatno + (forward ? -1 : 1) &&
8090 prevstate !=
nullptr &&
8094 jacscat[i].block(0, jmin, 4, jcnt) = jac.block(0, jmin, 4, jcnt);
8095 jacscat[i](4, 4) = jac(4, 4);
8097 calculateJac(jac, jacscat[i], jmin, jmax);
8101 Eigen::MatrixXd & derivmat = state->
derivatives();
8102 const int scatterPos = nperpars + 2 * i;
8104 derivmat.block<4, 2>(0, scatterPos) = (forward ? 1 : -1) * jacscat[i].block<4, 2>(0, 2);
8110 jcnt = jmax - jmin + 1;
8112 lp_bgn = forward ? nbremupstream : nbremupstream - 1;
8115 for (
int i = lp_bgn; forward ? (i < lp_end) : (i > lp_end); i += (forward ? 1 : -1)) {
8117 i == bremno + (forward ? -1 : 1) &&
8122 jacbrem[i].block(0, jmin, 4, jcnt) = jac.block(0, jmin, 4, jcnt);
8123 jacbrem[i](4, 4) = jac(4, 4);
8125 calculateJac(jac, jacbrem[i], jmin, jmax);
8129 Eigen::MatrixXd & derivmat = state->
derivatives();
8130 const int scatterPos = nperpars + 2 * nscats + i;
8132 derivmat.block<5, 1>(0, scatterPos) = (forward ? .001 : -.001) * jacbrem[i].block<5, 1>(0, 4);
8136 calculateJac(jac, jacvertex, 0, 4);
8140 Eigen::MatrixXd & derivmat = state->
derivatives();
8141 derivmat.block(0, 0, 4, nperpars) = jacvertex.block(0, 0, 4, nperpars);
8143 if (nperpars == 5) {
8144 derivmat.col(4).segment(0, 4) *= .001;
8145 derivmat(4, 4) = .001 * jacvertex(4, 4);
8151 (!trajectory.
prefit() || states[hitno]->materialEffects()->deltaE() == 0)
8153 scatno += (forward ? 1 : -1);
8157 states[hitno]->materialEffects() &&
8158 states[hitno]->materialEffects()->sigmaDeltaE() > 0
8160 bremno += (forward ? 1 : -1);
8163 prevstate = states[hitno].get();
8172 bool onlylocal)
const {
8178 std::vector<std::unique_ptr<GXFTrackState>> & states = trajectory.
trackStates();
8180 std::vector < int >
indices(states.size());
8182 int i = nstatesupstream;
8183 for (
int j = 0; j < (int) states.size(); j++) {
8184 if (j < nstatesupstream) {
8191 for (
int stateno = 0; stateno < (int) states.size(); stateno++) {
8192 if (stateno == 0 || stateno == nstatesupstream) {
8193 prevstate =
nullptr;
8196 std::unique_ptr<GXFTrackState> & state = states[
index];
8197 if (state->materialEffects() !=
nullptr) {
8198 prevstate = state.get();
8202 if (!state->hasTrackCovariance()) {
8203 state->zeroTrackCovariance();
8205 AmgMatrix(5, 5) & trackerrmat = state->trackCovariance();
8207 if ((prevstate !=
nullptr) &&
8211 Eigen::Matrix<double, 5, 5> & jac = state->jacobian();
8212 const AmgMatrix(5, 5)& prevcov = states[
indices[stateno - 1]]->trackCovariance();
8214 trackerrmat = jac * prevcov * jac.transpose();
8218 trackerrmat = derivatives * fullcovmat * derivatives.transpose();
8228 bool errorok =
true;
8229 for (
int i = 0; i < 5; i++) {
8232 && trackerrmat(i, i) > meascov(j, j)) {
8234 const double scale = std::sqrt(meascov(j, j) / trackerrmat(i, i));
8235 trackerrmat(i, i) = meascov(j, j);
8236 for (
int k = 0; k < 5; k++) {
8238 trackerrmat(k, i) *= scale;
8246 for (
int i = 0; i < 5; i++) {
8250 for (
int j = 0; j < 5; j++) {
8258 trackerrmat(4, 4) = 1e-20;
8262 state->trackParameters();
8264 std::optional<
AmgMatrix(5, 5)> trkerrmat;
8266 if (state->hasTrackCovariance()) {
8267 trkerrmat = (state->trackCovariance());
8269 trkerrmat = std::nullopt;
8272 const AmgVector(5) & tpars = tmptrackpar->parameters();
8273 std::unique_ptr<const TrackParameters> trackpar(
8279 std::move(trkerrmat))
8281 state->setTrackParameters(std::move(trackpar));
8284 if (errorok && trajectory.
nDOF() > 0) {
8285 fitQual =
m_updator->fullStateFitQuality(
8286 *state->trackParameters(),
8294 state->setFitQuality(fitQual);
8296 prevstate = state.get();
8300 std::optional<TransportJacobian>
8302 const EventContext& ctx,
8315 std::optional<TransportJacobian> jac = std::make_optional<TransportJacobian>(J);
8318 0.01, 0.01, 0.00001, 0.00001, 0.000000001
8329 for (
int i = 0; i < 5; i++) {
8332 if (thisdiscsurf && i == 1) {
8338 if (i == 0 && thiscylsurf) {
8339 vecminuseps[i] = -std::remainder(-vecminuseps[i], 2 *
M_PI * previousSurface.
bounds().
r());
8340 }
else if (i == 1 && thisdiscsurf) {
8341 vecpluseps[i] = -std::remainder(-vecpluseps[i], 2 *
M_PI);
8346 std::unique_ptr<const TrackParameters> parpluseps(
8356 const std::unique_ptr<const TrackParameters> parminuseps(
8367 std::unique_ptr<const TrackParameters> newparpluseps(
8378 std::unique_ptr<const TrackParameters> newparminuseps(
8393 if (newparpluseps ==
nullptr) {
8405 if (newparminuseps ==
nullptr) {
8417 if ((newparpluseps ==
nullptr) || (newparminuseps ==
nullptr)) {
8421 for (
int j = 0; j < 5; j++) {
8425 if (j == 0 && cylsurf) {
8427 }
else if (j == 1 && discsurf) {
8431 (*jac) (j, i) =
diff / (2 * eps[i]);
8444 (
"Configure the minimum number of Iterations via jobOptions");
8466 auto nmeas1 = pDataVector->size();
8467 const auto *pLastValue = (*pDataVector)[nmeas1 - 1];
8475 testrot =
static_cast<const RIO_OnTrack *
>(pLastValue);
8477 if (lastMeasIsCompetingRIO){
8479 testrot = &testcrot->rioOnTrack(0);
8483 if (testrot ==
nullptr) {
8484 const auto *pPenultimate = (*pDataVector)[nmeas1 - 2];
8487 if(penultimateIsRIO){
8488 testrot =
static_cast<const RIO_OnTrack *
>(pPenultimate);
8490 if (penultimateIsCompetingRIO){
8492 testrot = &testcrot->rioOnTrack(0);
8499 (testrot !=
nullptr) &&
8516 if (cond_obj ==
nullptr) {
8525 std::stringstream
msg;
8527 throw std::runtime_error(
msg.str());
8531 if (
cache.m_caloEntrance ==
nullptr) {
8534 if (geometry !=
nullptr) {
8535 cache.m_caloEntrance = geometry->trackingVolume(
"InDet::Containers::InnerDetector");
8543 if (
cache.m_caloEntrance ==
nullptr) {
8548 return cache.m_caloEntrance !=
nullptr;
8552 if (
cache.m_msEntrance ==
nullptr) {
8555 if (geometry !=
nullptr) {
8556 cache.m_msEntrance = geometry->trackingVolume(
"MuonSpectrometerEntrance");
8564 if (
cache.m_msEntrance ==
nullptr) {
8569 return cache.m_msEntrance !=
nullptr;
Scalar perp() const
perp method - perpendicular length
Scalar deltaPhi(const MatrixBase< Derived > &vec) const
Scalar mag() const
mag method
#define ATH_CHECK
Evaluate an expression and check for errors.
#define ATH_MSG_VERBOSE(x)
#define ATH_MSG_WARNING(x)
This class provides an interface to generate or decode an identifier for the upper levels of the dete...
double charge(const T &p)
std::vector< size_t > vec
#define AmgSymMatrix(dim)
#define AmgMatrix(rows, cols)
std::pair< std::vector< unsigned int >, bool > res
static const uint32_t nHits
void diff(const Jet &rJet1, const Jet &rJet2, std::map< std::string, double > varDiff)
Difference between jets - Non-Class function required by trigger.
size_t size() const
Number of registered mappings.
void getInitializedCache(MagField::AtlasFieldCache &cache) const
get B field cache for evaluation as a function of 2-d or 3-d position.
DataModel_detail::const_iterator< DataVector > const_iterator
Standard const_iterator.
const_iterator end() const noexcept
Return a const_iterator pointing past the end of the collection.
const_iterator begin() const noexcept
Return a const_iterator pointing at the beginning of the collection.
size_type size() const noexcept
Returns the number of elements in the collection.
bool empty() const noexcept
Returns true if the collection is empty.
bool is_valid() const
Check if id is in a valid state.
Class to hold geometrical description of a silicon detector element.
virtual bool isValid() override final
Can the handle be successfully dereferenced?
virtual std::span< T *const > arrayObjects()=0
Return all objects of the Array non-const we can still modify the T.
Base class for all CompetingRIOsOnTack implementations, extends the common MeasurementBase.
virtual const RIO_OnTrack & rioOnTrack(unsigned int) const =0
returns the RIO_OnTrack (also known as ROT) objects depending on the integer.
Bounds for a cylindrical Surface.
virtual double r() const override final
This method returns the radius.
double halflengthZ() const
This method returns the halflengthZ.
Class to describe a cylindrical detector layer for tracking, it inhertis from both,...
Class for a CylinderSurface in the ATLAS detector.
virtual const CylinderBounds & bounds() const override final
This method returns the CylinderBounds by reference (NoBounds is not possible for cylinder).
Bounds for a cylindrical Volume, the decomposeToSurfaces method creates a vector of up to 6 surfaces:
Class to describe the bounds for a planar DiscSurface.
double rMax() const
This method returns outer radius.
double rMin() const
This method returns inner radius.
Class to describe a disc-like detector layer for tracking, it inhertis from both, Layer base class an...
Class for a DiscSurface in the ATLAS detector.
const SurfaceBounds & bounds() const override final
This method returns the bounds by reference.
Access to distance solutions.
double second() const
Distance to second intersection solution along direction (for a cylinder surface).
int numberOfSolutions() const
Number of intersection solutions.
double first() const
Distance to first intersection solution along direction.
double sigmaDeltaE() const
returns the symmatric error
double deltaE() const
returns the
int numberDoF() const
returns the number of degrees of freedom of the overall track or vertex fit as integer
double chiSquared() const
returns the of the overall track fit
Status codes for track fitters.
@ ExtrapolationFailureDueToSmallMomentum
extrapolation failed due to small momentum
@ OutlierLogicFailure
outlier logic failed
@ ExtrapolationFailure
extrapolation failed
class that is similar to MaterialEffectsOnTrack, but has 'set' methods for more flexibility during tr...
void setSigmaDeltaE(double)
double deltaTheta() const
double sigmaDeltaEPos() const
double sigmaDeltaE() const
const Surface & associatedSurface() const
double sigmaDeltaPhi() const
void setScatteringSigmas(double, double)
double sigmaDeltaTheta() const
void setMaterialProperties(const MaterialProperties *)
Set the material properties of this material effects instance.
const MaterialProperties * materialProperties() const
bool getStateType(TrackStateOnSurface::TrackStateOnSurfaceType type) const
Retrieve the value of a specific type bit.
Amg::MatrixX & derivatives()
Eigen::Matrix< double, 5, 5 > & jacobian()
void setSinStereo(double)
const MeasurementBase * measurement(void)
TrackState::MeasurementType measurementType()
bool hasTrackCovariance(void) const
void setMeasurement(std::unique_ptr< const MeasurementBase >)
void setRecalibrated(bool)
void setMeasurementErrors(const double *)
const TrackParameters * trackParameters(void) const
double * measurementErrors()
GXFMaterialEffects * materialEffects()
void setTrackParameters(std::unique_ptr< const TrackParameters >)
const Surface & associatedSurface() const
Internal representation of the track, used in the track fit.
int numberOfOutliers() const
std::vector< double > & brems()
void setNumberOfScatterers(int)
int numberOfSiliconHits() const
std::vector< std::pair< const Layer *, const Layer * > > & upstreamMaterialLayers()
int numberOfFitParameters() const
int numberOfTRTPrecHits() const
Amg::MatrixX & weightedResidualDerivatives()
const std::vector< std::unique_ptr< GXFTrackState > > & trackStates() const
int numberOfTRTTubeHits() const
int numberOfTRTHits() const
int numberOfUpstreamBrems() const
double totalEnergyLoss() const
int numberOfUpstreamScatterers() const
int numberOfUpstreamStates() const
MagneticFieldProperties m_fieldprop
void setNumberOfPerigeeParameters(int)
const TrackParameters * referenceParameters()
void setNumberOfBrems(int)
std::vector< std::pair< double, double > > & scatteringSigmas()
std::pair< GXFTrackState *, GXFTrackState * > findFirstLastMeasurement(void)
void addMaterialState(std::unique_ptr< GXFTrackState >, int index=-1)
int numberOfScatterers() const
bool addMeasurementState(std::unique_ptr< GXFTrackState >, int index=-1)
std::vector< std::pair< double, double > > & scatteringAngles()
void addBasicState(std::unique_ptr< GXFTrackState >, int index=-1)
Amg::VectorX & residuals()
void setBrems(std::vector< double > &)
int numberOfPerigeeParameters() const
int numberOfBrems() const
void setOutlier(int, bool isoutlier=true)
void updateTRTHitCount(int index, float oldError)
void setScatteringAngles(std::vector< std::pair< double, double > > &)
void setReferenceParameters(std::unique_ptr< const TrackParameters >)
Gaudi::Property< int > m_maxoutliers
ToolHandle< IBoundaryCheckTool > m_boundaryCheckTool
std::unique_ptr< Track > makeTrack(const EventContext &ctx, Cache &, GXFTrajectory &, const ParticleHypothesis) const
void tryToConverge(const Cache &cache, GXFTrajectory &trajectory, const int it) const
void makeProtoState(Cache &, GXFTrajectory &, const TrackStateOnSurface *, int index=-1) const
Gaudi::Property< double > m_p
FitterStatusCode calculateTrackParameters(const EventContext &ctx, GXFTrajectory &, bool) const
Track * myfit(const EventContext &ctx, Cache &, GXFTrajectory &, const TrackParameters &, const RunOutlierRemoval runOutlier=false, const ParticleHypothesis matEffects=nonInteracting) const
void fillResidualsAndErrors(const EventContext &ctx, const Cache &cache, GXFTrajectory &trajectory, const int it, Amg::VectorX &b, int &bremno_maxbrempull, GXFTrackState *&state_maxbrempull) const
Gaudi::Property< double > m_scalefactor
Gaudi::Property< bool > m_rejectLargeNScat
std::variant< std::unique_ptr< const TrackParameters >, FitterStatusCode > updateEnergyLoss(const Surface &, const GXFMaterialEffects &, const TrackParameters &, double, int) const
ToolHandle< Trk::ITrkMaterialProviderTool > m_caloMaterialProvider
Gaudi::Property< bool > m_trtrecal
bool processTrkVolume(Cache &, const Trk::TrackingVolume *tvol) const
ToolHandle< IPropagator > m_propagator
Gaudi::Property< bool > m_straightlineprop
Gaudi::Property< bool > m_fiteloss
std::optional< GlobalChi2Fitter::TrackHoleCount > holeSearchProcess(const EventContext &ctx, const std::vector< std::reference_wrapper< GXFTrackState > > &states) const
Conduct a hole search between a list of states, possibly reusing existing information.
static std::optional< std::pair< Amg::Vector3D, double > > addMaterialFindIntersectionDisc(Cache &cache, const DiscSurface &surface, const TrackParameters ¶m1, const TrackParameters ¶m2, const ParticleHypothesis mat)
Find the intersection of a set of track parameters onto a disc surface.
ToolHandle< INavigator > m_navigator
void updateSystemWithMaxBremPull(GXFTrajectory &trajectory, const int bremno_maxbrempull, GXFTrackState *state_maxbrempull, Amg::SymMatrixX &a) const
SG::ReadHandleKey< Trk::ClusterSplitProbabilityContainer > m_clusterSplitProbContainer
virtual ~GlobalChi2Fitter()
ToolHandle< IRIO_OnTrackCreator > m_broadROTcreator
std::unique_ptr< GXFTrackState > makeTrackFindPerigee(const EventContext &, Cache &, GXFTrajectory &, const ParticleHypothesis) const
FitterStatusCode updateFitParameters(GXFTrajectory &, const Amg::VectorX &, const Amg::SymMatrixX &) const
Method to update peregee parameters, scattering angles, and brems.
void addMaterial(const EventContext &ctx, Cache &, GXFTrajectory &, const TrackParameters *, ParticleHypothesis) const
static void fillFirstLastMeasurement(Cache &cache, GXFTrajectory &trajectory)
PropagationResult calculateTrackParametersPropagateHelper(const EventContext &, const TrackParameters &, const GXFTrackState &, PropDirection, const MagneticFieldProperties &, bool, bool) const
Helper method that encapsulates calls to the propagator tool in the calculateTrackParameters() method...
ToolHandle< IUpdator > m_updator
static void fillBfromMeasurements(const Cache &cache, GXFTrajectory &trajectory, Amg::VectorX &b)
std::unique_ptr< const TrackParameters > makeTrackFindPerigeeParameters(const EventContext &, Cache &, GXFTrajectory &, const ParticleHypothesis) const
SG::ReadCondHandleKey< AtlasFieldCacheCondObj > m_field_cache_key
Gaudi::Property< bool > m_domeastrackpar
virtual int iterationsOfLastFit() const
Gaudi::Property< int > m_fixbrem
Track * backupCombinationStrategy(const EventContext &ctx, Cache &, const Track &, const Track &, GXFTrajectory &, std::vector< MaterialEffectsOnTrack > &) const
GlobalChi2Fitter(const std::string &, const std::string &, const IInterface *)
Track * mainCombinationStrategy(const EventContext &ctx, Cache &, const Track &, const Track &, GXFTrajectory &, std::vector< MaterialEffectsOnTrack > &) const
static void compensatePhiWeights(Cache &cache, GXFTrajectory &trajectory, Amg::SymMatrixX &a)
const TrackingGeometry * trackingGeometry(Cache &cache, const EventContext &ctx) const
Gaudi::Property< int > m_maxit
PropagationResult calculateTrackParametersPropagate(const EventContext &, const TrackParameters &, const GXFTrackState &, PropDirection, const MagneticFieldProperties &, bool, bool) const
Propagate onto a track state, collecting new track parameters, and optionally the Jacobian and possib...
ToolHandle< IExtrapolator > m_extrapolator
void makeProtoStateFromMeasurement(Cache &, GXFTrajectory &, const MeasurementBase *, const TrackParameters *trackpar=nullptr, bool isoutlier=false, int index=-1) const
std::optional< TransportJacobian > numericalDerivatives(const EventContext &ctx, const TrackParameters *, const Surface &, PropDirection, const MagneticFieldProperties &) const
void addIDMaterialFast(const EventContext &ctx, Cache &cache, GXFTrajectory &track, const TrackParameters *parameters, ParticleHypothesis part) const
A faster strategy for adding scatter material to tracks, works only for inner detector tracks.
static void fillAfromScatterers(GXFTrajectory &trajectory, Amg::SymMatrixX &a)
void updatePixelROTs(GXFTrajectory &, Amg::SymMatrixX &, Amg::VectorX &, const EventContext &evtctx) const
Update the Pixel ROT using the current trajectory/local track parameters.
Gaudi::Property< double > m_outlcut
Gaudi::Property< bool > m_calomat
Gaudi::Property< bool > m_fillderivmatrix
ToolHandle< IMaterialEffectsOnTrackProvider > m_calotoolparam
ToolHandle< IRIO_OnTrackCreator > m_ROTcreator
GXFTrajectory * runTrackCleanerSilicon(const EventContext &ctx, Cache &, GXFTrajectory &, Amg::SymMatrixX &, Amg::SymMatrixX &, Amg::VectorX &, bool) const
std::unique_ptr< const TrackParameters > makePerigee(Cache &, const TrackParameters &, const ParticleHypothesis) const
virtual std::unique_ptr< Track > fit(const EventContext &ctx, const PrepRawDataSet &, const TrackParameters &, const RunOutlierRemoval runOutlier=false, const ParticleHypothesis matEffects=nonInteracting) const override final
Gaudi::Property< double > m_chi2cut
void throwFailedToGetTrackingGeomtry() const
virtual Track * alignmentFit(const EventContext &ctx, AlignmentCache &, const Track &, const RunOutlierRemoval runOutlier=false, const ParticleHypothesis matEffects=Trk::nonInteracting) const override
Gaudi::Property< bool > m_extensioncuts
virtual void setMinIterations(int)
bool ensureValidEntranceCalo(const EventContext &ctx, Cache &cache) const
static void addMaterialGetLayers(Cache &cache, std::vector< std::pair< const Layer *, const Layer * > > &layers, std::vector< std::pair< const Layer *, const Layer * > > &uplayers, const std::vector< std::unique_ptr< GXFTrackState > > &states, GXFTrackState &first, GXFTrackState &last, const TrackParameters *refpar, bool hasmat)
Collect all possible layers that a given track could have passed through.
static std::optional< std::pair< Amg::Vector3D, double > > addMaterialFindIntersectionCyl(Cache &cache, const CylinderSurface &surface, const TrackParameters ¶m1, const TrackParameters ¶m2, const ParticleHypothesis mat)
Find the intersection of a set of track parameters onto a cylindrical surface.
void runTrackCleanerTRT(Cache &, GXFTrajectory &, Amg::SymMatrixX &, Amg::VectorX &, Amg::SymMatrixX &, bool, bool, int, const EventContext &ctx) const
static void calculateDerivatives(GXFTrajectory &)
Gaudi::Property< bool > m_decomposesegments
void calculateTrackErrors(GXFTrajectory &, Amg::SymMatrixX &, bool) const
static bool correctAngles(double &, double &)
bool ensureValidEntranceMuonSpectrometer(const EventContext &ctx, Cache &cache) const
bool isMuonTrack(const Track &) const
std::vector< std::unique_ptr< TrackParameters > > holesearchExtrapolation(const EventContext &ctx, const TrackParameters &src, const GXFTrackState &dst, PropDirection propdir) const
Helper method which performs an extrapolation with additional logic for hole search.
virtual StatusCode finalize() override
Gaudi::Property< bool > m_useCaloTG
std::vector< std::reference_wrapper< GXFTrackState > > holeSearchStates(GXFTrajectory &trajectory) const
Extracts a collection of track states which are important for hole search.
Track * fitIm(const EventContext &ctx, Cache &cache, const Track &inputTrack, const RunOutlierRemoval runOutlier, const ParticleHypothesis matEffects) const
Gaudi::Property< bool > m_redoderivs
void initFieldCache(const EventContext &ctx, Cache &cache) const
Initialize a field cache inside a fit cache object.
Gaudi::Property< int > m_maxitPixelROT
ToolHandle< IMaterialEffectsOnTrackProvider > m_calotool
ToolHandle< IMaterialEffectsUpdator > m_matupdator
void holeSearchHelper(const std::vector< std::unique_ptr< TrackParameters > > &hc, std::set< Identifier > &id_set, std::set< Identifier > &sct_set, TrackHoleCount &rv, bool count_holes, bool count_dead) const
Helper method for the hole search that does the actual counting of holes and dead modules.
ToolHandle< IResidualPullCalculator > m_residualPullCalculator
void addMaterialUpdateTrajectory(Cache &cache, GXFTrajectory &track, int offset, std::vector< std::pair< const Layer *, const Layer * > > &layers, const TrackParameters *ref1, const TrackParameters *ref2, ParticleHypothesis mat) const
Given layer information, probe those layers for scatterers and add them to a track.
ToolHandle< IEnergyLossUpdator > m_elosstool
void fillDerivatives(GXFTrajectory &traj) const
Gaudi::Property< bool > m_numderiv
SG::ReadCondHandleKey< TrackingGeometry > m_trackingGeometryReadKey
Gaudi::Property< bool > m_acceleration
static bool tryToWeightAfromMaterial(Cache &cache, GXFTrajectory &trajectory, Amg::SymMatrixX &a, const bool doDeriv, const int it, const double oldRedChi2, const double newRedChi2)
const AtlasDetectorID * m_DetID
static void fillAfromMeasurements(const Cache &cache, GXFTrajectory &trajectory, Amg::SymMatrixX &a)
virtual StatusCode initialize() override
static void makeTrackFillDerivativeMatrix(Cache &, GXFTrajectory &)
Gaudi::Property< bool > m_holeSearch
FitterStatusCode runIteration(const EventContext &ctx, Cache &cache, GXFTrajectory &trajectory, const int it, Amg::SymMatrixX &a, Amg::VectorX &b, Amg::SymMatrixX &lu, bool &doDeriv) const
Gaudi::Property< bool > m_createSummary
ToolHandle< IMultipleScatteringUpdator > m_scattool
Interface class IPropagators It inherits from IAlgTool.
LayerIndex for the identification of layers in a simplified detector geometry of Cylinders and Discs.
int value() const
layerIndex expressed in an integer
double oppositePostFactor() const
Return method for post update material description of the Layer along normalvector.
double alongPreFactor() const
Return method for pre update material description of the Layer along normalvector.
double alongPostFactor() const
Return method for post update material description of the Layer along normalvector.
double oppositePreFactor() const
Return method for pre update material description of the Layer along normalvector.
Base Class for a Detector Layer in the Tracking realm.
virtual const Surface & surfaceRepresentation() const =0
Transforms the layer into a Surface representation for extrapolation.
const MaterialProperties * fullUpdateMaterialProperties(const TrackParameters &par) const
getting the MaterialProperties back - for full update
const LayerMaterialProperties * layerMaterialProperties() const
getting the LayerMaterialProperties including full/pre/post update
int parameterKey() const
Identifier key for matrix expansion/reduction.
bool contains(ParamDefs par) const
The simple check for the clients whether the parameter is contained.
magnetic field properties to steer the behavior of the extrapolation
base class to integrate material effects on Trk::Track in a flexible way.
const Surface & associatedSurface() const
returns the surface to which these m.eff. are associated.
double thicknessInX0() const
returns the actually traversed material .
@ MATERIAL_EFFECTS_ON_TRACK
virtual MaterialEffectsDerivedType derivedType() const =0
Returns the concrete derived type.
represents the full description of deflection and e-loss of a track in material.
const EnergyLoss * energyLoss() const
returns the energy loss object.
const ScatteringAngles * scatteringAngles() const
returns the MCS-angles object.
Material with information about thickness of material.
float thicknessInX0() const
Return the radiationlength fraction.
float thickness() const
Return the thickness in mm.
This class is the pure abstract base class for all fittable tracking measurements.
virtual MeasurementBase * clone() const =0
Pseudo-Constructor.
const LocalParameters & localParameters() const
Interface method to get the LocalParameters.
virtual const Surface & associatedSurface() const =0
Interface method to get the associated Surface.
virtual bool type(MeasurementBaseType::Type type) const =0
Interface method checking the type.
virtual const Amg::Vector3D & globalPosition() const =0
Interface method to get the global Position.
const Amg::MatrixX & localCovariance() const
Interface method to get the localError.
const Amg::Vector3D & momentum() const
Access method for the momentum.
virtual ParametersBase< DIM, T > * clone() const override=0
clone method for polymorphic deep copy
const Amg::Vector3D & position() const
Access method for the position.
virtual const Surface & associatedSurface() const override=0
Access to the Surface associated to the Parameters.
double pT() const
Access method for transverse momentum.
Class describing the Line to which the Perigee refers to.
virtual DistanceSolution straightLineDistanceEstimate(const Amg::Vector3D &pos, const Amg::Vector3D &dir) const override final
fast straight line distance evaluation to Surface
virtual Surface::ChargedTrackParametersUniquePtr createUniqueTrackParameters(double l1, double l2, double phi, double theta, double qop, std::optional< AmgSymMatrix(5)> cov=std::nullopt) const override final
Use the Surface as a ParametersBase constructor, from local parameters - charged.
Class for a planaer rectangular or trapezoidal surface in the ATLAS detector.
const Amg::Vector2D & localPosition() const
return the local position reference
virtual bool type(PrepRawDataType type) const
Interface method checking the type.
const Amg::MatrixX & localCovariance() const
return const ref to the error matrix
Class to handle RIO On Tracks ROT) for InDet and Muons, it inherits from the common MeasurementBase.
virtual const Surface & associatedSurface() const override=0
returns the surface for the local to global transformation
virtual const Trk::PrepRawData * prepRawData() const =0
returns the PrepRawData (also known as RIO) object to which this RIO_OnTrack is associated.
Identifier identify() const
return the identifier -extends MeasurementBase
virtual const Amg::Vector3D & globalPosition() const override=0
Interface method to get the global Position.
@ Biased
RP with track state including the hit.
represents a deflection of the track caused through multiple scattering in material.
Base class for all TrackSegment implementations, extends the common MeasurementBase.
const MeasurementBase * measurement(unsigned int) const
returns the Trk::MeasurementBase objects depending on the integer
unsigned int numberOfMeasurementBases() const
Return the number of contained Trk::MeasurementBase (s).
Class for a StraightLineSurface in the ATLAS detector to describe dirft tube and straw like detectors...
virtual double r() const =0
Interface method for the maximal extension or the radius.
virtual BoundsType type() const =0
Return the bounds type - for persistency optimization.
Abstract Base Class for tracking surfaces.
virtual DistanceSolution straightLineDistanceEstimate(const Amg::Vector3D &pos, const Amg::Vector3D &dir) const =0
fast straight line distance evaluation to Surface
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.
virtual const Amg::Vector3D & normal() const
Returns the normal vector of the Surface (i.e.
const Amg::Transform3D & transform() const
Returns HepGeom::Transform3D by reference.
Identifier associatedDetectorElementIdentifier() const
return Identifier of the associated Detector Element
virtual constexpr SurfaceType type() const =0
Returns the Surface type to avoid dynamic casts.
virtual const SurfaceBounds & bounds() const =0
Surface Bounds method.
const Trk::Layer * associatedLayer() const
return the associated Layer
const Amg::Vector3D & center() const
Returns the center position of the Surface.
Contains information about the 'fitter' of this track.
bool trackProperties(const TrackProperties &property) const
Access methods for track properties.
@ GlobalChi2Fitter
Track's from Thijs' global chi^2 fitter.
@ Unknown
Track fitter not defined.
ParticleHypothesis particleHypothesis() const
Returns the particle hypothesis used for Track fitting.
const TrackFitter & trackFitter() const
Access methods for track fitter.
@ BremFit
A brem fit was performed on this track.
@ StraightTrack
A straight track.
@ BremFitSuccessful
A brem fit was performed on this track and this fit was successful.
@ SlimmedTrack
A slimmed track.
represents the track state (measurement, material, fit parameters and quality) at a surface.
const MeasurementBase * measurementOnTrack() const
returns MeasurementBase const overload
const TrackParameters * trackParameters() const
return ptr to trackparameters const overload
bool type(const TrackStateOnSurfaceType type) const
Use this method to find out if the TSoS is of a certain type: i.e.
@ Measurement
This is a measurement, and will at least contain a Trk::MeasurementBase.
@ Perigee
This represents a perigee, and so will contain a Perigee object only.
@ BremPoint
This represents a brem point on the track, and so will contain TrackParameters and MaterialEffectsBas...
@ Outlier
This TSoS contains an outlier, that is, it contains a MeasurementBase/RIO_OnTrack which was not used ...
@ InertMaterial
This represents inert material, and so will contain MaterialEffectsBase.
@ Scatterer
This represents a scattering point on the track, and so will contain TrackParameters and MaterialEffe...
@ CaloDeposit
This TSOS contains a CaloEnergy object.
const MaterialEffectsBase * materialEffectsOnTrack() const
return material effects const overload
const Trk::Surface & surface() const
return associated surface
const Trk::TrackStates * trackStateOnSurfaces() const
return a pointer to a const DataVector of const TrackStateOnSurfaces.
const DataVector< const MeasurementBase > * measurementsOnTrack() const
return a pointer to a vector of MeasurementBase (NOT including any that come from outliers).
const DataVector< const TrackParameters > * trackParameters() const
Return a pointer to a vector of TrackParameters.
const TrackInfo & info() const
Returns a const ref to info of a const tracks.
const Perigee * perigeeParameters() const
return Perigee.
const FitQuality * fitQuality() const
return a pointer to the fit quality const-overload
The TrackingGeometry class is the owner of the constructed TrackingVolumes.
Full Volume description used in Tracking, it inherits from Volume to get the geometrical structure,...
const LayerArray * confinedLayers() const
Return the subLayer array.
LayerIntersection< Amg::Vector3D > closestMaterialLayer(const Amg::Vector3D &gp, const Amg::Vector3D &dir, PropDirection pDir=alongMomentum, const BoundaryCheck &bchk=true) const
Return the closest layer with material description.
const TrackingVolumeArray * confinedVolumes() const
Return the subLayer array.
std::vector< std::shared_ptr< BoundarySurface< TrackingVolume > > > & boundarySurfaces()
Method to return the BoundarySurfaces.
This is the base class for all tracking detector elements with read-out relevant information.
virtual Identifier identify() const =0
Identifier.
double chi2(TH1 *h0, TH1 *h1)
std::string head(std::string s, const std::string &pattern)
head of a string
Eigen::Matrix< double, 3, 3 > RotationMatrix3D
Eigen::Matrix< double, Eigen::Dynamic, Eigen::Dynamic > MatrixX
Dynamic Matrix - dynamic allocation.
Eigen::Matrix< double, Eigen::Dynamic, Eigen::Dynamic > SymMatrixX
Eigen::Affine3d Transform3D
Eigen::Matrix< double, 3, 1 > Vector3D
Eigen::Matrix< double, Eigen::Dynamic, 1 > VectorX
Dynamic Vector - dynamic allocation.
double getDistance(const xAOD::Vertex *vtx1, const xAOD::Vertex *vtx2)
@ PseudoMeasurementOnTrack
constexpr double mass[PARTICLEHYPOTHESES]
the array of masses
MeasurementType
enum describing the flavour of MeasurementBase
Ensure that the ATLAS eigen extensions are properly loaded.
PropDirection
PropDirection, enum for direction of the propagation.
std::vector< const MeasurementBase * > MeasurementSet
vector of fittable measurements
ParametersT< TrackParametersDim, Charged, PerigeeSurface > Perigee
ParametersT< TrackParametersDim, Charged, StraightLineSurface > AtaStraightLine
@ DeadElement
outside the element
bool consistentSurfaces(U)
std::unique_ptr< T > unique_clone(const T *v)
bool RunOutlierRemoval
switch to toggle quality processing after fit
@ NoField
Field is set to 0., 0., 0.,.
@ FullField
Field is set to be realistic, but within a given Volume.
@ z
global position (cartesian)
@ u
Enums for curvilinear frames.
@ loc2
generic first and second local coordinate
std::pair< double, ParamDefs > DefinedParameter
Typedef to of a std::pair<double, ParamDefs> to identify a passed-through double as a specific type o...
ParticleHypothesis
Enumeration for Particle hypothesis respecting the interaction with material.
std::pair< long int, long int > indices
BinnedArray< TrackingVolume > TrackingVolumeArray
simply for the eye
const Amg::Vector3D & position() const
Method to retrieve the position of the Intersection.
ParametersBase< TrackParametersDim, Charged > TrackParameters
ParametersT< TrackParametersDim, Charged, PlaneSurface > AtaPlane
const IIntersectionCache * cache() const
Retrieve the associated cache block, if it exists.
std::vector< const PrepRawData * > PrepRawDataSet
vector of clusters and drift circles
AmgSymMatrix(5) &GXFTrackState
@ numberOfSCTHoles
number of Holes in both sides of a SCT module
@ numberOfSCTDeadSensors
number of TRT hits
@ numberOfPixelHoles
number of pixels which have a ganged ambiguity.
@ numberOfPixelDeadSensors
number of pixel hits with broad errors (width/sqrt(12))
void sort(typename DataModel_detail::iterator< DVL > beg, typename DataModel_detail::iterator< DVL > end)
Specialization of sort for DataVector/List.
void stable_sort(DataModel_detail::iterator< DVL > beg, DataModel_detail::iterator< DVL > end)
Specialization of stable_sort for DataVector/List.
void reverse(typename DataModel_detail::iterator< DVL > beg, typename DataModel_detail::iterator< DVL > end)
Specialization of reverse for DataVector/List.
double deltaPhi(double phiA, double phiB)
delta Phi in range [-pi,pi[
static void objVectorDeleter(const std::vector< const T * > *ptr)
static constexpr std::array< ParamDefs, 6 > pardef
Constructor.