56#include "CLHEP/Units/SystemOfUnits.h"
70#include <Eigen/StdVector>
81 return distsol.
first();
84 std::abs(distsol.
first()) < std::abs(distsol.
second()) ?
93 bool correctAngles(
double&
phi,
double&
theta) {
110 double limitInversePValue(
double qOverP) {
111 const double magnitude = std::abs(qOverP);
113 constexpr double maxP{100.*10e6*
MeV};
114 constexpr double minP{1.e-3*
MeV};
115 constexpr double lo {1./maxP};
116 constexpr double hi {1./minP};
117 const double limited = std::clamp(magnitude, lo, hi);
118 return std::copysign(limited, qOverP);
121 std::pair<const Trk::TrackParameters *, const Trk::TrackParameters *> getFirstLastIdPar(
const Trk::Track & track) {
127 while ((firstidpar ==
nullptr) && parit !=
track.trackParameters()->end()) {
129 ((**parit).covariance() !=
nullptr) &&
138 parit =
track.trackParameters()->end();
142 ((**parit).covariance() !=
nullptr) &&
147 }
while ((lastidpar ==
nullptr) && parit !=
track.trackParameters()->begin());
149 return std::make_pair(firstidpar, lastidpar);
164 std::abs(
a.parameters()[0] -
b.parameters()[0]) < e &&
165 std::abs(
a.parameters()[1] -
b.parameters()[1]) < e &&
166 std::abs(
a.parameters()[2] -
b.parameters()[2]) < e
177 calculateJac(Eigen::Matrix<double, 5, 5> &jac,
178 Eigen::Matrix<double, 5, 5> &out,
179 int jmin,
int jmax) {
183 out.block(0, 0, 4, jmin).setZero();
187 out.block(0, jmax + 1, 4, 5 - (jmax + 1)).setZero();
190 out(4, 4) = jac(4, 4);
200 std::pair<double, double> principalComponentAnalysis2x2(
const Amg::MatrixX & mat) {
201 const double trace =
mat(0, 0) +
mat(1, 1);
202 const double diagonalProduct =
mat(0, 0) *
mat(1, 1);
203 const double mat01Sq =
mat(0, 1) *
mat(0, 1);
204 const double discriminant = std::sqrt(trace * trace - 4. * (diagonalProduct - mat01Sq));
206 const double eigenValueSmall = 0.5 * (trace -
discriminant);
207 const double stereoAngle = 0.5 * std::asin(2 *
mat(0, 1) / (-discriminant));
209 return std::make_pair(eigenValueSmall, stereoAngle);
215 const std::string & t,
216 const std::string & n,
244 ATH_MSG_ERROR(
"Hole search requested but no boundary check tool provided.");
245 return StatusCode::FAILURE;
270 ATH_MSG_WARNING(
"FillDerivativeMatrix option selected, switching off acceleration!");
294 ATH_MSG_ERROR(
"Hole search requested but track summaries are disabled.");
295 return StatusCode::FAILURE;
300 return StatusCode::SUCCESS;
306 if (m_fit_status[
S_FITS] > 0) {
309 <<
" track fits failed because of a matrix inversion failure");
311 <<
" tracks were rejected by the outlier logic");
313 <<
" track fits failed because of a propagation failure");
315 <<
" track fits failed because of an invalid angle (theta/phi)");
317 <<
" track fits failed because the fit did not converge");
319 <<
" tracks did not pass the chi^2 cut");
321 <<
" tracks were killed by the energy loss update");
324 return StatusCode::SUCCESS;
329 std::unique_ptr<Track>
331 const EventContext& ctx,
337 ATH_MSG_DEBUG(
"--> entering GlobalChi2Fitter::fit(Track,Track,)");
345 !
cache.m_field_cache.solenoidOn() && !
cache.m_field_cache.toroidOn()
352 const Track *indettrack = firstismuon ? &intrk2 : &intrk1;
353 const Track *muontrack = firstismuon ? &intrk1 : &intrk2;
355 bool measphi =
false;
362 rot = &crot->rioOnTrack(0);
372 const double dotprod2 = measdir.dot(
375 if (std::abs(dotprod1) < .5 && std::abs(dotprod2) < .5) {
383 auto [firstidpar, lastidpar] = getFirstLastIdPar(*indettrack);
385 if ((firstidpar ==
nullptr) || (lastidpar ==
nullptr)) {
389 std::unique_ptr<const TrackParameters> parforcalo =
unique_clone(firstismuon ? firstidpar : lastidpar);
391 if (!
cache.m_field_cache.solenoidOn()) {
392 const AmgVector(5) & newpars = parforcalo->parameters();
394 parforcalo=parforcalo->associatedSurface().createUniqueTrackParameters(
395 newpars[0], newpars[1], newpars[2], newpars[3], 1 / 5000., std::nullopt
399 std::vector < MaterialEffectsOnTrack > calomeots;
402 calomeots =
m_calotool->extrapolationSurfacesAndEffects(
406 parforcalo->associatedSurface(),
419 if (calomeots.empty()) {
424 std::unique_ptr<Track> track;
427 cache.m_calomat =
false;
428 const bool tmp2 =
cache.m_extmat;
429 const bool tmp4 =
cache.m_idmat;
434 const double qoverpid = measperid !=
nullptr ? measperid->parameters()[
Trk::qOverP] : 0;
435 const double qoverpmuon = measpermuon !=
nullptr ? measpermuon->parameters()[
Trk::qOverP] : 0;
437 const AmgSymMatrix(5) * errmatid = measperid !=
nullptr ? measperid->covariance() :
nullptr;
438 const AmgSymMatrix(5) * errmatmuon = measpermuon !=
nullptr ? measpermuon->covariance() :
nullptr;
442 (errmatid !=
nullptr) &&
443 (errmatmuon !=
nullptr) &&
448 const double piderror = std::sqrt((*errmatid) (4, 4)) / (qoverpid * qoverpid);
449 const double pmuonerror = std::sqrt((*errmatmuon) (4, 4)) / (qoverpmuon * qoverpmuon);
450 const double energyerror = std::sqrt(
451 calomeots[1].energyLoss()->sigmaDeltaE() *
452 calomeots[1].energyLoss()->sigmaDeltaE() + piderror * piderror +
453 pmuonerror * pmuonerror
457 (std::abs(calomeots[1].energyLoss()->deltaE()) -
458 std::abs(1 / qoverpid) + std::abs(1 / qoverpmuon)
461 ATH_MSG_DEBUG(
"Changing from measured to parametrized energy loss");
466 parforcalo->associatedSurface(),
471 if (calomeots.empty()) {
479 bool firstfitwasattempted =
false;
482 if (!caloEntranceIsValid) {
487 (!
cache.m_field_cache.toroidOn() && !
cache.m_field_cache.solenoidOn()) ||
489 cache.m_getmaterialfromtrack &&
493 qoverpid * qoverpmuon > 0
498 if (
cache.m_fit_status[
S_FITS] == (
unsigned int) (nfits + 1)) {
499 firstfitwasattempted =
true;
504 (track ==
nullptr) &&
505 !firstfitwasattempted &&
506 (
cache.m_field_cache.toroidOn() ||
cache.m_field_cache.solenoidOn())
512 trajectory = trajectory2;
516 bool pseudoupdated =
false;
518 if (track !=
nullptr) {
519 for (std::unique_ptr<GXFTrackState> & pseudostate : trajectory.
trackStates()) {
520 if (pseudostate ==
nullptr) {
531 if ((pseudostate ==
nullptr) || pseudostate->fitQuality().chiSquared() < 10) {
536 const std::unique_ptr<const TrackParameters> updpar(
m_updator->removeFromState(
538 pseudostate->measurement()->localParameters(),
539 pseudostate->measurement()->localCovariance()
542 if (updpar ==
nullptr) {
547 covMatrix(0, 0) = 100;
549 std::unique_ptr<const PseudoMeasurementOnTrack> newpseudo = std::make_unique<const PseudoMeasurementOnTrack>(
553 std::move(covMatrix),
557 pseudostate->setMeasurement(std::move(newpseudo));
559 errors[0] = errors[2] = errors[3] = errors[4] = -1;
561 pseudostate->setMeasurementErrors(errors);
562 pseudoupdated =
true;
567 cache.m_matfilled =
true;
573 *track->perigeeParameters(),
578 cache.m_matfilled =
false;
584 if (track !=
nullptr) {
585 track->info().addPatternReco(intrk1.
info());
586 track->info().addPatternReco(intrk2.
info());
590 cache.m_calomat = tmp;
591 cache.m_extmat = tmp2;
592 cache.m_idmat = tmp4;
597 const EventContext& ctx,
599 const Track & intrk1,
600 const Track & intrk2,
602 std::vector<MaterialEffectsOnTrack> & calomeots
604 ATH_MSG_DEBUG(
"--> entering GlobalChi2Fitter::mainCombinationStrategy");
609 const Track *indettrack = firstismuon ? &intrk2 : &intrk1;
610 const Track *muontrack = firstismuon ? &intrk1 : &intrk2;
612 auto [tmpfirstidpar, tmplastidpar] = getFirstLastIdPar(*indettrack);
613 std::unique_ptr<const TrackParameters> firstidpar =
unique_clone(tmpfirstidpar);
614 std::unique_ptr<const TrackParameters> lastidpar =
unique_clone(tmplastidpar);
616 if ((firstidpar ==
nullptr) || (lastidpar ==
nullptr)) {
630 std::unique_ptr<const TrackParameters> tp_closestmuon =
nullptr;
632 while (closestmuonmeas ==
nullptr) {
633 closestmuonmeas =
nullptr;
636 if ((**tsosit).measurementOnTrack() !=
nullptr) {
637 closestmuonmeas = (**tsosit).measurementOnTrack();
639 if (thispar !=
nullptr) {
640 const AmgVector(5) & parvec = thispar->parameters();
642 parvec[0], parvec[1], parvec[2], parvec[3], parvec[4], std::nullopt
656 std::unique_ptr<const TrackParameters> tmppar;
659 if ((tp_closestmuon !=
nullptr) && msEntranceIsValid) {
664 std::unique_ptr<const std::vector<const TrackStateOnSurface *>> matvec;
666 if (tmppar !=
nullptr) {
667 const Surface & associatedSurface = tmppar->associatedSurface();
668 std::unique_ptr<Surface> muonsurf =
nullptr;
674 const double radius = cylbounds->
r();
676 muonsurf = std::make_unique<CylinderSurface>(trans, radius + 1, hlength);
680 const double newz = (
681 associatedSurface.
center().
z() > 0 ?
682 associatedSurface.
center().z() + 1 :
683 associatedSurface.
center().z() - 1
687 associatedSurface.
center().x(),
688 associatedSurface.
center().y(),
692 trans.translation() << newpos;
695 const double rmin = discbounds->
rMin();
696 const double rmax = discbounds->
rMax();
697 muonsurf = std::make_unique<DiscSurface>(trans, rmin, rmax);
701 if (muonsurf !=
nullptr) {
713 std::vector<const TrackStateOnSurface *> tmp_matvec;
715 if ((matvec !=
nullptr) && !matvec->empty()) {
716 tmp_matvec = *matvec;
717 delete tmp_matvec.back();
718 tmp_matvec.pop_back();
720 for (
auto & i : tmp_matvec) {
739 if (tmppar ==
nullptr) {
753 if (tmppar ==
nullptr) {
757 AmgVector(5) newpars = tmppar->parameters();
762 const double oldp = std::abs(1 / newpars[
Trk::qOverP]);
763 const double newp2 = oldp * oldp + (!firstismuon ? 2 : -2) * de * std::sqrt(mass * mass + oldp * oldp) + de * de;
770 tp_closestmuon=tmppar->associatedSurface().createUniqueTrackParameters(
771 newpars[0], newpars[1], newpars[2], newpars[3], newpars[4], std::nullopt
787 for (; itStates != endState; ++itStates) {
792 const bool tmpgetmat =
cache.m_getmaterialfromtrack;
794 if ((*itStates)->materialEffectsOnTrack() !=
nullptr) {
796 cache.m_extmat =
false;
798 cache.m_idmat =
false;
801 const auto *
const pBaseMEOT = (*itStates)->materialEffectsOnTrack();
805 const auto *
const pMEOT =
static_cast<const MaterialEffectsOnTrack *
>((*itStates)->materialEffectsOnTrack());
806 if ((pMEOT->scatteringAngles() ==
nullptr) or (pMEOT->energyLoss() ==
nullptr)) {
807 cache.m_getmaterialfromtrack =
true;
813 cache.m_getmaterialfromtrack = tmpgetmat;
820 trajectory.
trackStates().back()->setTrackParameters(
nullptr);
823 std::unique_ptr<const TrackParameters> firstscatpar;
824 std::unique_ptr<const TrackParameters> lastscatpar;
827 double newqoverpid = 0;
830 const double de = std::abs(calomeots[1].energyLoss()->deltaE());
831 const double sigmade = std::abs(calomeots[1].energyLoss()->sigmaDeltaE());
833 const double pbefore = std::abs(1 / firstidpar->parameters()[
Trk::qOverP]);
834 const double pafter = std::abs(1 / tp_closestmuon->parameters()[
Trk::qOverP]);
835 const double elosspull = (pbefore - pafter - de) / sigmade;
837 if (std::abs(elosspull) > 10) {
838 if (elosspull > 10) {
839 newqoverpid = 1 / (de + pafter + 10 * sigmade);
841 newqoverpid = 1 / (de + pafter - 10 * sigmade);
844 if (tp_closestmuon->parameters()[
Trk::qOverP] * newqoverpid < 0) {
848 const AmgVector(5) & newpar = firstidpar->parameters();
849 firstidpar=firstidpar->associatedSurface().createUniqueTrackParameters(
850 newpar[0], newpar[1], newpar[2], newpar[3], newqoverpid, std::nullopt
858 if (lastidpar ==
nullptr) {
864 *(firstismuon ? tp_closestmuon.get() : lastidpar.get()),
865 calomeots[0].associatedSurface(),
872 if (firstscatpar ==
nullptr) {
878 *(firstismuon ? firstidpar : tp_closestmuon),
879 calomeots[2].associatedSurface(),
886 if (lastscatpar ==
nullptr) {
890 std::optional<TransportJacobian> jac1;
891 std::optional<TransportJacobian> jac2;
892 std::unique_ptr<const TrackParameters> elosspar;
894 double firstscatphi = 0;
895 double secondscatphi = 0;
896 double firstscattheta = 0;
897 double secondscattheta = 0;
898 double muonscatphi = 0;
899 double muonscattheta = 0;
901 const TrackParameters *idscatpar = !firstismuon ? firstscatpar.get() : lastscatpar.get();
902 const TrackParameters *muonscatpar = firstismuon ? firstscatpar.get() : lastscatpar.get();
904 newqoverpid = idscatpar->parameters()[
Trk::qOverP];
906 const Amg::Vector3D calosegment = lastscatpar->position() - firstscatpar->position();
909 muonscattheta = calosegment.theta() - muonscatpar->parameters()[
Trk::theta];
912 for (
int i = 0; i < 2; i++) {
913 std::unique_ptr<const TrackParameters> tmpelosspar;
914 AmgVector(5) params1 = muonscatpar->parameters();
923 params1[0], params1[1], params1[2], params1[3], params1[4], std::nullopt
944 calomeots[1].associatedSurface(),
950 if ((tmpelosspar ==
nullptr) || (jac1 == std::nullopt)) {
959 const AmgVector(5) & newpars = tmpelosspar->parameters();
960 const std::unique_ptr<const TrackParameters> elosspar2(tmpelosspar->associatedSurface().createUniqueTrackParameters(
961 newpars[0], newpars[1], newpars[2], newpars[3], newqoverpid, std::nullopt
965 elosspar = std::move(tmpelosspar);
968 std::unique_ptr<const TrackParameters> scat2(
m_propagator->propagateParameters(
972 calomeots[0].associatedSurface() :
973 calomeots[2].associatedSurface(),
986 calomeots[0].associatedSurface() :
987 calomeots[2].associatedSurface(),
995 if ((scat2 ==
nullptr) || (jac2 == std::nullopt)) {
1000 for (
int j = 0; j < 5; j++) {
1001 for (
int k = 0; k < 5; k++) {
1003 for (
int l = 0; l < 5; l++) {
1004 jac3[j][k] += (*jac2) (j, l) * (*jac1) (l, k);
1013 jac4(0, 0) = jac3[0][2];
1014 jac4(1, 1) = jac3[1][3];
1015 jac4(0, 1) = jac3[0][3];
1016 jac4(1, 0) = jac3[1][2];
1018 jac4 = jac4.inverse();
1030 discsurf =
static_cast<const Trk::DiscSurface *
>(&scat2->associatedSurface());
1032 if (cylsurf !=
nullptr) {
1033 dloc1 = -std::remainder(-dloc1, 2 *
M_PI * cylsurf->
bounds().
r());
1036 if (discsurf !=
nullptr) {
1037 dloc2 = -std::remainder(-dloc2, 2 *
M_PI);
1040 double dphi = jac4(0, 0) * dloc1 + jac4(0, 1) * dloc2;
1041 double dtheta = jac4(1, 0) * dloc1 + jac4(1, 1) * dloc2;
1047 muonscatphi += dphi;
1048 muonscattheta += dtheta;
1051 const double idscattheta = idscatpar->parameters()[
Trk::theta] - (scat2->parameters()[
Trk::theta] + dtheta);
1054 firstscatphi = muonscatphi;
1055 secondscatphi = idscatphi;
1056 firstscattheta = muonscattheta;
1057 secondscattheta = idscattheta;
1059 firstscatphi = -idscatphi;
1060 secondscatphi = -muonscatphi;
1061 firstscattheta = -idscattheta;
1062 secondscattheta = -muonscattheta;
1065 if (i == 1 &&
cache.m_field_cache.toroidOn() && !firstismuon) {
1066 AmgVector(5) params2 = scat2->parameters();
1074 firstscatpar=scat2->associatedSurface().createUniqueTrackParameters(
1075 params2[0], params2[1], params2[2], params2[3], params2[4], std::nullopt
1077 idscatpar = firstscatpar.get();
1081 *
cache.m_caloEntrance,
1085 if (startPar !=
nullptr) {
1094 rot.col(2) = trackdir;
1097 trans.linear().matrix() << rot;
1098 trans.translation() << startPar->position() - .1 * trackdir;
1109 if (curvlinpar !=
nullptr) {
1110 startPar.reset(curvlinpar);
1114 firstscatpar = std::move(scat2);
1118 std::unique_ptr<GXFMaterialEffects> firstscatmeff = std::make_unique<GXFMaterialEffects>(calomeots[0]);
1119 std::unique_ptr<GXFMaterialEffects> elossmeff = std::make_unique<GXFMaterialEffects>(calomeots[1]);
1120 std::unique_ptr<GXFMaterialEffects> secondscatmeff = std::make_unique<GXFMaterialEffects>(calomeots[2]);
1122 const double pull1 = std::abs(firstscatphi / firstscatmeff->sigmaDeltaPhi());
1123 const double pull2 = std::abs(secondscatphi / secondscatmeff->sigmaDeltaPhi());
1126 for (
auto & i : tmp_matvec) {
1131 firstscatmeff->setScatteringAngles(firstscatphi, firstscattheta);
1132 secondscatmeff->setScatteringAngles(secondscatphi, secondscattheta);
1135 elossmeff->setdelta_p(1000 * (lastscatpar->parameters()[
Trk::qOverP] - newqoverpid));
1137 elossmeff->setdelta_p(1000 * (newqoverpid - firstscatpar->parameters()[
Trk::qOverP]));
1140 elossmeff->setSigmaDeltaE(calomeots[1].energyLoss()->sigmaDeltaE());
1142 trajectory.
addMaterialState(std::make_unique<GXFTrackState>(std::move(firstscatmeff), std::move(firstscatpar)), -1);
1143 trajectory.
addMaterialState(std::make_unique<GXFTrackState>(std::move(elossmeff), std::move(elosspar)), -1);
1144 trajectory.
addMaterialState(std::make_unique<GXFTrackState>(std::move(secondscatmeff), std::move(lastscatpar)), -1);
1147 for (
auto & i : tmp_matvec) {
1154 if (startPar ==
nullptr) {
1159 (pull1 > 5 || pull2 > 5) &&
1165 bool largegap =
false;
1166 double previousz = 0;
1168 for (itStates2 = beginStates2; itStates2 != endState2; ++itStates2) {
1173 if (meff !=
nullptr) {
1186 cache.m_extmat =
false;
1188 cache.m_idmat =
false;
1194 !(itStates2 == beginStates2 || itStates2 == beginStates2 + 1) &&
1203 itStates2 == endState2 - 1 &&
1206 std::abs(tpar->
position().z()) < 13000
1208 std::unique_ptr<const TrackParameters> pseudopar(tpar->
clone());
1210 covMatrix(0, 0) = 100;
1212 std::unique_ptr<const PseudoMeasurementOnTrack> newpseudo = std::make_unique<const PseudoMeasurementOnTrack>(
1214 std::move(covMatrix),
1215 pseudopar->associatedSurface()
1218 std::unique_ptr<GXFTrackState> pseudostate = std::make_unique<GXFTrackState>(std::move(newpseudo), std::move(pseudopar));
1222 errors[0] = errors[2] = errors[3] = errors[4] = -1;
1225 pseudostate->setMeasurementErrors(errors);
1232 std::abs(trajectory.
trackStates().back()->position().z()) > 20000 &&
1233 std::abs(previousz) < 12000
1239 previousz = trajectory.
trackStates().back()->position().z();
1257 const EventContext& ctx,
1259 const Track & intrk1,
1260 const Track & intrk2,
1262 std::vector<MaterialEffectsOnTrack> & calomeots
1264 ATH_MSG_DEBUG(
"--> entering GlobalChi2Fitter::backupCombinationStrategy");
1267 const Track *indettrack = firstismuon ? &intrk2 : &intrk1;
1280 if (pParametersVector->size() > 1)
1281 firstidpar = (*pParametersVector)[1];
1283 firstidpar = pParametersVector->back();
1285 std::unique_ptr<const TrackParameters> lastidpar =
nullptr;
1286 if ((firstidpar !=
nullptr) && (
cache.m_caloEntrance !=
nullptr))
1290 if (lastidpar ==
nullptr) {
1291 lastidpar.reset(pParametersVector->back()->clone());
1294 std::unique_ptr < const TrackParameters > firstscatpar;
1295 std::unique_ptr < const TrackParameters > lastscatpar;
1296 std::unique_ptr < const TrackParameters > elosspar;
1301 lastidpar->position(),
1302 lastidpar->momentum(),
1304 lastidpar->position()
1311 calomeots[0].associatedSurface(),
1317 if (!firstscatpar) {
1321 const std::unique_ptr<const TrackParameters> tmppar(
1325 calomeots[1].associatedSurface(),
1338 const double oldp = std::abs(1 / tmppar->parameters()[
Trk::qOverP]);
1339 const double de = std::abs(calomeots[1].energyLoss()->deltaE());
1341 double newp2 = oldp * oldp - 2 * de * std::sqrt(mass * mass + oldp * oldp) + de * de;
1347 const double newqoverp =
sign / std::sqrt(newp2);
1349 const AmgVector(5) & pars = tmppar->parameters();
1352 tmppar->associatedSurface().createUniqueTrackParameters(
1353 pars[0], pars[1], pars[2], pars[3], newqoverp, std::nullopt
1359 calomeots[2].associatedSurface(),
1372 calomeots[2].associatedSurface(),
1386 calomeots[1].associatedSurface(),
1397 const double sign = (elosspar->parameters()[
Trk::qOverP] < 0) ? -1 : 1;
1398 const double newqoverp =
sign /
1399 (1. / std::abs(elosspar->parameters()[
Trk::qOverP]) +
1400 std::abs(calomeots[1].energyLoss()->deltaE()));
1402 const AmgVector(5) & pars = elosspar->parameters();
1404 std::unique_ptr<const TrackParameters>
const tmppar(
1405 elosspar->associatedSurface().createUniqueTrackParameters(
1406 pars[0], pars[1], pars[2], pars[3], newqoverp, std::nullopt
1413 calomeots[0].associatedSurface(),
1420 if (!firstscatpar) {
1425 for (; itStates != endState; ++itStates) {
1433 if ((*itStates)->materialEffectsOnTrack() !=
nullptr) {
1435 cache.m_idmat =
false;
1446 std::unique_ptr<GXFMaterialEffects> firstscatmeff = std::make_unique<GXFMaterialEffects>(calomeots[0]);
1447 std::unique_ptr<GXFMaterialEffects> elossmeff = std::make_unique<GXFMaterialEffects>(calomeots[1]);
1448 std::unique_ptr<GXFMaterialEffects> secondscatmeff = std::make_unique<GXFMaterialEffects>(calomeots[2]);
1452 sigmadp = calomeots[1].energyLoss()->sigmaDeltaE();
1453 elossmeff->setSigmaDeltaE(sigmadp);
1456 elossmeff->setdelta_p(dp);
1458 trajectory.
addMaterialState(std::make_unique<GXFTrackState>(std::move(firstscatmeff), std::move(firstscatpar)), -1);
1459 trajectory.
addMaterialState(std::make_unique<GXFTrackState>(std::move(elossmeff), std::move(elosspar)), -1);
1460 trajectory.
addMaterialState(std::make_unique<GXFTrackState>(std::move(secondscatmeff), std::move(lastscatpar)), -1);
1463 const Surface *triggersurf1 =
nullptr;
1464 const Surface *triggersurf2 =
nullptr;
1468 bool seenmdt =
false;
1469 bool mdtbetweenphihits =
false;
1473 itStates2 = (!firstismuon ? beginStates2 : endState - 1);
1474 itStates2 != (!firstismuon ? endState2 : beginStates - 1);
1475 (!firstismuon ? ++itStates2 : --itStates2)
1478 ((*itStates2)->measurementOnTrack() ==
nullptr) ||
1483 const auto *
const pMeasurement = (*itStates2)->measurementOnTrack();
1484 const Surface *surf = &pMeasurement->associatedSurface();
1488 if (isCompetingRIOsOnTrack) {
1490 rot = &crot->rioOnTrack(0);
1493 rot =
static_cast<const RIO_OnTrack *
>(pMeasurement);
1496 if ((rot !=
nullptr) &&
m_DetID->is_mdt(rot->
identify()) && (triggersurf1 !=
nullptr)) {
1500 (rot !=
nullptr) && (
1507 const double dotprod1 = measdir.dot(
Amg::Vector3D(0, 0, 1));
1510 const bool measphi = std::abs(dotprod1) <= .5 && std::abs(dotprod2) <= .5;
1514 (*itStates2)->trackParameters() !=
nullptr ?
1515 (*itStates2)->trackParameters()->position() :
1517 if (triggersurf1 !=
nullptr) {
1518 triggerpos2 = thispos;
1519 triggersurf2 = surf;
1521 mdtbetweenphihits =
true;
1524 triggerpos1 = thispos;
1525 triggersurf1 = surf;
1531 double mdttrig1 = 999999;
1532 double mdttrig2 = 999999;
1533 const Surface *mdtsurf1 =
nullptr;
1534 const Surface *mdtsurf2 =
nullptr;
1537 itStates2 = (!firstismuon ? beginStates2 : endState - 1);
1538 itStates2 != (!firstismuon ? endState2 : beginStates - 1);
1539 (!firstismuon ? ++itStates2 : --itStates2)
1541 const Surface *surf =
nullptr;
1543 ((*itStates2)->measurementOnTrack() !=
nullptr) &&
1546 surf = &(*itStates2)->measurementOnTrack()->associatedSurface();
1549 if (surf ==
nullptr) {
1552 const auto *
const pThisMeasurement = (*itStates2)->measurementOnTrack();
1556 if (isCompetingRioOnTrack) {
1558 rot = &crot->rioOnTrack(0);
1561 rot =
static_cast<const RIO_OnTrack *
>(pThisMeasurement);
1567 (*itStates2)->trackParameters() !=
nullptr ?
1568 (*itStates2)->trackParameters()->position() :
1569 pThisMeasurement->globalPosition();
1570 if (triggerpos1.mag() > 1 && (globpos - triggerpos1).mag() < mdttrig1) {
1571 mdttrig1 = (globpos - triggerpos1).
mag();
1574 if (triggerpos2.mag() > 1 && (globpos - triggerpos2).mag() < mdttrig2) {
1575 mdttrig2 = (globpos - triggerpos2).
mag();
1582 std::vector<GXFTrackState *> outlierstates;
1583 std::vector<GXFTrackState *> outlierstates2;
1585 outlierstates.reserve(10);
1586 outlierstates2.reserve(10);
1588 std::unique_ptr<PseudoMeasurementOnTrack> newpseudo;
1590 for (itStates2 = beginStates2; itStates2 != endState2; ++itStates2) {
1591 const auto *
const pMeasurement{(*itStates2)->measurementOnTrack()};
1593 const bool isStraightLine =
1594 pMeasurement !=
nullptr ?
1601 (newpseudo ==
nullptr) && (
1602 (itStates2 == beginStates2 || itStates2 == beginStates2 + 1) &&
1603 std::abs(pMeasurement->globalPosition().z()) < 10000
1606 std::unique_ptr<const TrackParameters> par2;
1607 if (((*itStates2)->trackParameters() !=
nullptr) && nphi > 99) {
1608 par2.reset((*itStates2)->trackParameters()->clone());
1620 if (par2 ==
nullptr) {
1624 covMatrix(0, 0) = 100;
1626 newpseudo = std::make_unique<PseudoMeasurementOnTrack>(
1628 std::move(covMatrix),
1629 par2->associatedSurface()
1632 std::unique_ptr<GXFTrackState> firstpseudo = std::make_unique<GXFTrackState>(std::move(newpseudo), std::move(par2));
1636 errors[0] = errors[2] = errors[3] = errors[4] = -1;
1639 firstpseudo->setMeasurementErrors(errors);
1640 firstpseudostate = firstpseudo.get();
1646 if (isPseudo && !firstismuon) {
1650 if ((**itStates2).materialEffectsOnTrack() !=
nullptr) {
1652 cache.m_idmat =
false;
1660 ((**itStates2).measurementOnTrack() !=
nullptr) &&
1661 &(**itStates2).measurementOnTrack()->associatedSurface() == triggersurf1 &&
1662 (mdtsurf1 !=
nullptr)
1664 std::unique_ptr<Amg::Transform3D> transf = std::make_unique<Amg::Transform3D>(mdtsurf1->
transform());
1666 transf->translation() << triggerpos1;
1669 covMatrix(0, 0) = 100;
1671 std::unique_ptr<const PseudoMeasurementOnTrack> newpseudo = std::make_unique<const PseudoMeasurementOnTrack>(
1675 std::unique_ptr<GXFTrackState> pseudostate1 = std::make_unique<GXFTrackState>(std::move(newpseudo),
nullptr);
1679 errors[0] = errors[2] = errors[3] = errors[4] = -1;
1682 pseudostate1->setMeasurementErrors(errors);
1683 outlierstates2.push_back(pseudostate1.get());
1688 ((**itStates2).measurementOnTrack() !=
nullptr) &&
1689 &(**itStates2).measurementOnTrack()->associatedSurface() == triggersurf2 &&
1690 mdtbetweenphihits &&
1691 (mdtsurf2 !=
nullptr)
1693 std::unique_ptr<Amg::Transform3D> transf = std::make_unique<Amg::Transform3D>(mdtsurf2->
transform());
1694 transf->translation() << triggerpos2;
1697 covMatrix(0, 0) = 100;
1699 std::unique_ptr<const PseudoMeasurementOnTrack> newpseudo = std::make_unique<const PseudoMeasurementOnTrack>(
1703 std::unique_ptr<GXFTrackState> pseudostate2 = std::make_unique<GXFTrackState>(std::move(newpseudo),
nullptr);
1707 errors[0] = errors[2] = errors[3] = errors[4] = -1;
1710 pseudostate2->setMeasurementErrors(errors);
1712 outlierstates2.push_back(pseudostate2.get());
1728 outlierstates.push_back(trajectory.
trackStates().back().get());
1736 Track *track =
nullptr;
1740 cache.m_matfilled =
true;
1741 const bool tmpacc =
cache.m_acceleration;
1742 cache.m_acceleration =
false;
1744 const std::unique_ptr<Trk::Track> tmp_track(
1746 cache.m_acceleration = tmpacc;
1748 cache.m_matfilled =
false;
1752 std::abs(trajectory.
residuals().tail<1>()(0) / trajectory.
errors().tail<1>()(0)) > 10
1758 if (firstpseudostate !=
nullptr) {
1761 covMatrix(0, 0) = 100;
1763 std::unique_ptr<const PseudoMeasurementOnTrack> newpseudo = std::make_unique<const PseudoMeasurementOnTrack>(
1765 std::move(covMatrix),
1773 for (
const auto & i : outlierstates2) {
1779 for (
const auto & i : outlierstates) {
1787 itStates = (firstismuon ? beginStates2 : endState - 1);
1788 itStates != (firstismuon ? endState2 : beginStates - 1);
1789 (firstismuon ? ++itStates : --itStates)
1802 cache.m_matfilled =
false;
1808 std::unique_ptr<Track>
1810 const EventContext& ctx,
1811 const Track& inputTrack,
1815 ATH_MSG_DEBUG(
"--> entering GlobalChi2Fitter::fit(Track,)");
1828 return std::unique_ptr<Track>(
1829 fitIm(ctx,
cache, inputTrack, runOutlier, matEffects));
1834 AlignmentCache& alignCache,
1835 const Track &inputTrack,
1842 alignCache.m_derivMatrix.reset();
1843 alignCache.m_fullCovarianceMatrix.reset();
1844 alignCache.m_iterationsOfLastFit = 0;
1847 fitIm(ctx,
cache, inputTrack, runOutlier, matEffects);
1848 if(newTrack !=
nullptr){
1849 if(
cache.m_derivmat.size() != 0)
1850 alignCache.m_derivMatrix = std::make_unique<Amg::MatrixX>(
cache.m_derivmat);
1851 if(
cache.m_fullcovmat.size() != 0)
1852 alignCache.m_fullCovarianceMatrix = std::make_unique<Amg::MatrixX>(
cache.m_fullcovmat);
1853 alignCache.m_iterationsOfLastFit =
cache.m_lastiter;
1860 const EventContext& ctx,
1862 const Track& inputTrack,
1867 ATH_MSG_DEBUG(
"--> entering GlobalChi2Fitter::fit(Track,,)");
1883 ATH_MSG_WARNING(
"Track without track parameters, cannot perform fit");
1891 if (minpar ==
nullptr) {
1895 const bool tmpgetmat =
cache.m_getmaterialfromtrack;
1901 cache.m_getmaterialfromtrack =
false;
1909 const Surface *firsthitsurf =
nullptr;
1910 const Surface *lasthitsurf =
nullptr;
1912 bool hasmuon =
false;
1913 bool iscombined =
false;
1914 bool seenphimeas =
false;
1918 for (; itStates != endState; ++itStates) {
1919 const auto *
const pMeasurement = (**itStates).measurementOnTrack();
1921 (pMeasurement ==
nullptr) &&
1922 ((**itStates).materialEffectsOnTrack() !=
nullptr) &&
1928 if (pMeasurement !=
nullptr) {
1929 const Surface *surf = &pMeasurement->associatedSurface();
1938 if (firsthitsurf ==
nullptr) {
1939 firsthitsurf = surf;
1942 if (
m_DetID->is_indet(hitid)) {
1947 if ((**itStates).trackParameters() !=
nullptr) {
1948 lastidpar = (**itStates).trackParameters();
1949 if (firstidpar ==
nullptr) {
1950 firstidpar = lastidpar;
1956 const double dotprod1 = measdir.dot(
Amg::Vector3D(0, 0, 1));
1958 if (std::abs(dotprod1) < .5 && std::abs(dotprod2) < .5) {
1960 if (std::abs(surf->
center().z()) > 13000) {
1963 if (surf->
center().perp() > 9000 && std::abs(surf->
center().z()) < 13000) {
1975 if (iscombined && seenphimeas && (phiem || phibo)) {
1986 (hasmuon ||
cache.m_acceleration)
1988 cache.m_matfilled =
true;
1991 if (firstidpar == lastidpar) {
1992 firstidpar = lastidpar =
nullptr;
1997 !
cache.m_matfilled &&
2002 (firstidpar !=
nullptr)
2011 const bool tmpacc =
cache.m_acceleration;
2013 const bool tmpsirecal =
cache.m_sirecal;
2014 std::unique_ptr<Track> tmptrack =
nullptr;
2018 cache.m_fiteloss =
true;
2019 cache.m_sirecal =
false;
2022 cache.m_asymeloss =
true;
2025 tmptrack.reset(
myfit(ctx,
cache, trajectory, *minpar,
false, matEffects));
2026 cache.m_sirecal = tmpsirecal;
2028 if (tmptrack ==
nullptr) {
2029 cache.m_matfilled =
false;
2030 cache.m_getmaterialfromtrack = tmpgetmat;
2031 cache.m_acceleration = tmpacc;
2032 cache.m_fiteloss = tmpfiteloss;
2037 bool isbrem =
false;
2039 unsigned int n_brem=0;
2041 for (std::unique_ptr<GXFTrackState> & state : trajectory.
trackStates()) {
2044 if (meff !=
nullptr) {
2050 const double p = 1. / std::abs(layerpars->parameters()[
Trk::qOverP] - .0005 * meff->
delta_p());
2052 std::optional<Amg::Vector2D> locpos(state->associatedSurface().globalToLocal(layerpars->
position()));
2053 const Amg::Vector3D layerNormal(state->associatedSurface().normal(*locpos));
2054 double costracksurf = 1.;
2056 costracksurf = std::abs(layerNormal.dot(layerpars->
momentum().unit()));
2058 const double oldde = meff->
deltaE();
2060 std::unique_ptr<EnergyLoss> eloss;
2061 double sigmascat = 0;
2063 if (matprop !=
nullptr) {
2064 eloss = std::make_unique<EnergyLoss>(
2065 m_elosstool->energyLoss(*matprop, p, 1. / costracksurf,
2067 sigmascat = std::sqrt(
m_scattool->sigmaSquare(*matprop, p, 1. / costracksurf, matEffects));
2069 if (eloss !=
nullptr) {
2074 sigmascat = std::sqrt(
m_scattool->sigmaSquare(tmpprop, p, 1. / costracksurf, matEffects));
2078 sigmascat / std::sin(layerpars->parameters()[
Trk::theta]),
2092 }
else if (eloss !=
nullptr) {
2103 refpars[0], refpars[1], refpars[2], refpars[3], refpars[4], std::nullopt
2107 cache.m_matfilled =
true;
2111 trajectory.
brems().clear();
2113 trajectory.
brems().resize(1);
2114 trajectory.
brems()[0] = bremdp;
2117 cache.m_asymeloss =
false;
2125 std::unique_ptr<Track> track(
myfit(ctx,
cache, trajectory, *minpar, runOutlier, matEffects));
2127 bool pseudoupdated =
false;
2129 if ((track !=
nullptr) && hasid && hasmuon) {
2130 for (std::unique_ptr<GXFTrackState> & pseudostate : trajectory.
trackStates()) {
2132 (pseudostate ==
nullptr) ||
2134 pseudostate->fitQuality().chiSquared() < 10
2140 const std::unique_ptr<const TrackParameters> updpar(
m_updator->removeFromState(
2142 pseudostate->measurement()->localParameters(),
2143 pseudostate->measurement()->localCovariance()
2146 if (updpar ==
nullptr) {
2151 covMatrix(0, 0) = 100;
2153 std::unique_ptr<const PseudoMeasurementOnTrack> newpseudo = std::make_unique<const PseudoMeasurementOnTrack>(
2155 std::move(covMatrix),
2159 pseudostate->setMeasurement(std::move(newpseudo));
2161 errors[0] = errors[2] = errors[3] = errors[4] = -1;
2163 pseudostate->setMeasurementErrors(errors);
2164 pseudoupdated =
true;
2167 if (pseudoupdated) {
2169 cache.m_matfilled =
true;
2170 track.reset(
myfit(ctx,
cache, trajectory, *track->perigeeParameters(),
false,
muon));
2171 cache.m_matfilled =
false;
2175 cache.m_matfilled =
false;
2176 cache.m_getmaterialfromtrack = tmpgetmat;
2177 cache.m_acceleration = tmpacc;
2178 cache.m_fiteloss = tmpfiteloss;
2180 if (track !=
nullptr) {
2183 track->info().addPatternReco(old_info);
2186 return track.release();
2189 std::unique_ptr<Track>
2191 const EventContext& ctx,
2197 ATH_MSG_DEBUG(
"--> entering GlobalChi2Fitter::fit(PRDS,TP,)");
2200 for (
const auto *prd : prds) {
2201 const Surface & prdsurf = (*prd).detectorElement()->surface((*prd).identify());
2206 plsurf =
static_cast < const PlaneSurface *
>(&prdsurf);
2213 if ((slsurf ==
nullptr) && (plsurf ==
nullptr)) {
2214 ATH_MSG_ERROR(
"Surface is neither PlaneSurface nor StraightLineSurface!");
2219 }
else if (slsurf !=
nullptr) {
2228 }
else if (plsurf !=
nullptr) {
2229 if (param.covariance() !=
nullptr) {
2251 if (rot !=
nullptr) {
2252 rots.push_back(rot);
2256 std::unique_ptr<Track> track =
2257 fit(ctx, rots, param, runOutlier, matEffects);
2259 for (
const auto *rot : rots) {
2266 std::unique_ptr<Track>
2268 const EventContext& ctx,
2269 const Track& inputTrack,
2274 ATH_MSG_DEBUG(
"--> entering GlobalChi2Fitter::fit(Track,Meas'BaseSet,,)");
2289 if (minpar ==
nullptr) {
2299 const bool old_reintoutl =
cache.m_reintoutl;
2300 cache.m_reintoutl =
false;
2301 const bool tmpasymeloss =
cache.m_asymeloss;
2304 cache.m_asymeloss =
true;
2307 for (; itStates != endState; ++itStates) {
2312 (trajectory.
trackStates().back()->materialEffects() !=
nullptr) &&
2313 trajectory.
trackStates().back()->materialEffects()->sigmaDeltaE() > 50.001
2315 trajectory.
trackStates().back()->materialEffects()->setKink(
true);
2319 cache.m_reintoutl = old_reintoutl;
2321 for (
const auto & measBase : addMeasColl) {
2322 if (measBase ==
nullptr) {
2323 ATH_MSG_WARNING(
"There is an empty MeasurementBase object in the track! Skip this object..");
2331 std::unique_ptr<Track> track(
myfit(ctx,
cache, trajectory, *minpar, runOutlier, matEffects));
2332 cache.m_asymeloss = tmpasymeloss;
2334 if (track !=
nullptr) {
2335 const double oldqual =
2340 const double newqual =
2341 track->fitQuality()->numberDoF() != 0 ?
2342 track->fitQuality()->chiSquared() / track->fitQuality()->numberDoF() :
2345 if (
m_extensioncuts && runOutlier && newqual > 2 && newqual > 2 * oldqual) {
2349 track.reset(
nullptr);
2353 if (track !=
nullptr) {
2362 std::unique_ptr<Track>
2364 const EventContext& ctx,
2370 ATH_MSG_DEBUG(
"--> entering GlobalChi2Fitter::fit(Track,PRDS,)");
2374 for (
const auto *prd : prds) {
2375 const Surface & prdsurf = (*prd).detectorElement()->surface((*prd).identify());
2378 std::unique_ptr<const TrackParameters>
const trackparForCorrect(
2394 rot =
m_ROTcreator->correct(*prd, *trackparForCorrect, ctx);
2397 if (rot !=
nullptr) {
2398 rots.push_back(rot);
2402 std::unique_ptr<Track> track =
fit(ctx,intrk, rots, runOutlier, matEffects);
2404 for (
const auto *rot : rots) {
2412 const EventContext& ctx,
2418 ATH_MSG_DEBUG(
"--> entering GlobalChi2Fitter::fit(Meas'BaseSet,,)");
2431 for (
const auto *itSet : rots) {
2432 if (itSet ==
nullptr) {
2433 ATH_MSG_WARNING(
"There is an empty MeasurementBase object in the track! Skip this object..");
2439 std::unique_ptr<const TrackParameters> startpar(param.
clone());
2442 matEffects ==
muon &&
2445 cache.m_matfilled =
true;
2448 myfit(ctx,
cache, trajectory, *startpar, runOutlier, matEffects);
2450 cache.m_matfilled =
false;
2464 covMatrix(0, 0) = 100;
2466 std::unique_ptr<const PseudoMeasurementOnTrack> newpseudo = std::make_unique<const PseudoMeasurementOnTrack>(
2468 std::move(covMatrix),
2472 trajectory.
trackStates().front()->setMeasurement(std::move(newpseudo));
2477 covMatrix(0, 0) = 100;
2479 std::unique_ptr<const PseudoMeasurementOnTrack> newpseudo = std::make_unique<const PseudoMeasurementOnTrack>(
2481 std::move(covMatrix),
2485 trajectory.
trackStates().back()->setMeasurement(std::move(newpseudo));
2492 refpars[0], refpars[1], refpars[2], refpars[3], refpars[4], std::nullopt
2497 myfit(ctx,
cache, trajectory, *startpar, runOutlier, matEffects);
2499 cache.m_matfilled =
true;
2508 refpars[0], refpars[1], refpars[2], refpars[3], refpars[4], std::nullopt
2515 firstpar = trajectory.
trackStates().front()->trackParameters();
2518 covMatrix(0, 0) = 100;
2520 std::unique_ptr<const PseudoMeasurementOnTrack> newpseudo = std::make_unique<const PseudoMeasurementOnTrack>(
2522 std::move(covMatrix),
2526 trajectory.
trackStates().front()->setMeasurement(std::move(newpseudo));
2528 errors[0] = errors[2] = errors[3] = errors[4] = -1;
2530 trajectory.
trackStates().front()->setMeasurementErrors(errors);
2534 lastpar = trajectory.
trackStates().back()->trackParameters();
2537 covMatrix(0, 0) = 100;
2539 std::unique_ptr<const PseudoMeasurementOnTrack> newpseudo = std::make_unique<const PseudoMeasurementOnTrack>(
2541 std::move(covMatrix),
2545 trajectory.
trackStates().back()->setMeasurement(std::move(newpseudo));
2547 errors[0] = errors[2] = errors[3] = errors[4] = -1;
2549 trajectory.
trackStates().back()->setMeasurementErrors(errors);
2553 std::unique_ptr<Track> track;
2555 if (startpar !=
nullptr) {
2556 track.reset(
myfit(ctx,
cache, trajectory, *startpar, runOutlier, matEffects));
2559 if (track !=
nullptr) {
2562 cache.m_matfilled =
false;
2579 ) &&
cache.m_getmaterialfromtrack
2589 std::unique_ptr<GXFMaterialEffects> newmeff;
2597 newmeff = std::make_unique<GXFMaterialEffects>(*meff);
2601 const double sigmascat = std::sqrt(
m_scattool->sigmaSquare(
2622 newmeff = std::make_unique<GXFMaterialEffects>(newmeot);
2638 ((**trajectory.
trackStates().rbegin()).trackParameters() !=
nullptr)
2640 const double delta_p = 1000 * (
2642 (**trajectory.
trackStates().rbegin()).trackParameters()->
2646 newmeff->setdelta_p(delta_p);
2682 seg =
static_cast<const Segment *
>(measbase);
2692 for (
int i = 0; i <
imax; i++) {
2695 std::unique_ptr<GXFTrackState> ptsos = std::make_unique<GXFTrackState>(
2696 std::unique_ptr<const MeasurementBase>(measbase2->
clone()),
2697 std::unique_ptr<const TrackParameters>(newtrackpar !=
nullptr ? newtrackpar->
clone() :
nullptr)
2700 double sinstereo = 0;
2702 errors[0] = errors[1] = errors[2] = errors[3] = errors[4] = -1;
2713 bool measphi =
false;
2716 bool rotated =
false;
2719 if (
m_DetID->is_pixel(hitid)) {
2721 }
else if (
m_DetID->is_sct(hitid)) {
2722 if (covmat.cols() != 1 && covmat(1, 0) != 0) {
2726 }
else if (
m_DetID->is_trt(hitid)) {
2736 }
else if (
m_DetID->is_mdt(hitid)) {
2738 }
else if (
m_DetID->is_tgc(hitid)) {
2743 }
else if (
m_DetID->is_csc(hitid)) {
2745 }
else if (
m_DetID->is_mm(hitid)) {
2747 }
else if (
m_DetID->is_stgc(hitid)) {
2753 const auto [covEigenValueSmall, covStereoAngle] = principalComponentAnalysis2x2(covmat);
2754 errors[0] = std::sqrt(covEigenValueSmall);
2755 sinstereo = std::sin(covStereoAngle);
2757 errors[0] = std::sqrt(covmat(0, 0));
2759 errors[1] = std::sqrt(covmat(1, 1));
2770 const double dotprod1 = measdir.dot(
Amg::Vector3D(0, 0, 1));
2772 if (std::abs(dotprod1) < .5 && std::abs(dotprod2) < .5) {
2780 int param_index = 0;
2782 errors[0] = std::sqrt(covmat(0, 0));
2787 errors[1] = std::sqrt(covmat(param_index, param_index));
2792 errors[2] = std::sqrt(covmat(param_index, param_index));
2797 errors[3] = std::sqrt(covmat(param_index, param_index));
2802 errors[4] = std::sqrt(covmat(param_index, param_index));
2823 ptsos->setMeasurementErrors(errors);
2824 ptsos->setSinStereo(sinstereo);
2825 ptsos->setMeasurementType(hittype);
2826 ptsos->setMeasuresPhi(measphi);
2828 if (isoutlier && !
cache.m_reintoutl) {
2847 if (tvol ==
nullptr) {
2851 const Trk::BinnedArray < Trk::Layer > *confinedLayers = tvol->
confinedLayers();
2854 if (confinedLayers !=
nullptr) {
2856 for (
const auto & layer : confinedLayers->
arrayObjects()) {
2858 if (layer !=
nullptr) {
2863 if ((layIndex.
value() == 0) || (layer->layerMaterialProperties() ==
nullptr)) {
2874 disclay =
static_cast<const DiscLayer *
>(layer);
2877 if (disclay !=
nullptr) {
2878 if (disclay->
center().z() < 0) {
2879 cache.m_negdiscs.push_back(disclay);
2881 cache.m_posdiscs.push_back(disclay);
2883 }
else if (cyllay !=
nullptr) {
2884 cache.m_barrelcylinders.push_back(cyllay);
2894 for (
size_t ib = 0 ; ib < bsurf.size(); ++ib) {
2895 const Layer *layer = bsurf[ib]->surfaceRepresentation().materialLayer();
2897 if (layer ==
nullptr)
continue;
2901 if ((layIndex.
value() == 0) || (layer->layerMaterialProperties() ==
nullptr)) {
2908 cylsurf =
static_cast<const CylinderSurface *
>(&layer->surfaceRepresentation());
2913 discsurf =
static_cast<const DiscSurface *
>(&layer->surfaceRepresentation());
2915 if (discsurf !=
nullptr) {
2917 discsurf->
center().z() < 0 &&
2918 std::find(
cache.m_negdiscs.begin(),
cache.m_negdiscs.end(), layer) ==
cache.m_negdiscs.end()
2920 cache.m_negdiscs.push_back(layer);
2922 discsurf->
center().z() > 0 &&
2923 std::find(
cache.m_posdiscs.begin(),
cache.m_posdiscs.end(), layer) ==
cache.m_posdiscs.end()
2925 cache.m_posdiscs.push_back(layer);
2928 (cylsurf !=
nullptr) &&
2929 std::find(
cache.m_barrelcylinders.begin(),
cache.m_barrelcylinders.end(), layer) ==
cache.m_barrelcylinders.end()
2931 cache.m_barrelcylinders.push_back(layer);
2934 if ((cylsurf ==
nullptr) && (discsurf ==
nullptr)) {
2941 if (confinedVolumes !=
nullptr) {
2942 for (
const auto & volume : confinedVolumes->
arrayObjects()) {
2943 if (volume !=
nullptr) {
2959 std::optional<std::pair<Amg::Vector3D, double>>
2972 const double * pos = parforextrap.
position().data();
2974 cache.m_field_cache.getFieldZR(pos, field);
2975 const double sinphi = std::sin(parforextrap.parameters()[
Trk::phi0]);
2976 const double cosphi = std::cos(parforextrap.parameters()[
Trk::phi0]);
2977 const double sintheta = std::sin(parforextrap.parameters()[
Trk::theta]);
2978 const double costheta = std::cos(parforextrap.parameters()[
Trk::theta]);
2982 const double r = (std::abs(currentqoverp) > 1e-10) ? -sintheta / (currentqoverp * 300. * field[2]) : 1e6;
2983 const double xc = parforextrap.
position().x() -
r * sinphi;
2984 const double yc = parforextrap.
position().y() +
r * cosphi;
2985 const double phi0 = std::atan2(parforextrap.
position().y() - yc, parforextrap.
position().x() - xc);
2986 const double z0 = parforextrap.
position().z();
2987 const double delta_s = (surf.
center().
z() -
z0) / costheta;
2988 const double delta_phi = delta_s * sintheta /
r;
2989 const double x = xc + std::abs(
r) * std::cos(
phi0 + delta_phi);
2990 const double y = yc + std::abs(
r) * std::sin(
phi0 + delta_phi);
2992 const double perp = intersect.perp();
2999 const double costracksurf = std::abs(costheta);
3001 return std::make_pair(intersect, costracksurf);
3004 std::optional<std::pair<Amg::Vector3D, double>>
3019 const double * pos = parforextrap.
position().data();
3021 cache.m_field_cache.getFieldZR(pos, field);
3022 const double sinphi = std::sin(parforextrap.parameters()[
Trk::phi0]);
3023 const double cosphi = std::cos(parforextrap.parameters()[
Trk::phi0]);
3024 const double sintheta = std::sin(parforextrap.parameters()[
Trk::theta]);
3025 const double costheta = std::cos(parforextrap.parameters()[
Trk::theta]);
3026 const double tantheta = std::tan(parforextrap.parameters()[
Trk::theta]);
3027 const double r = (std::abs(currentqoverp) > 1e-10) ? -sintheta / (currentqoverp * 300. * field[2]) : 1e6;
3028 const double xc = parforextrap.
position().x() -
r * sinphi;
3029 const double yc = parforextrap.
position().y() +
r * cosphi;
3030 const double phi0 = std::atan2(parforextrap.
position().y() - yc, parforextrap.
position().x() - xc);
3031 const double z0 = parforextrap.
position().z();
3032 const double d = xc * xc + yc * yc;
3033 const double rcyl = surf.
bounds().
r();
3034 double mysqrt = ((
r + rcyl) * (
r + rcyl) - d) * (d - (
r - rcyl) * (
r - rcyl));
3040 mysqrt = std::sqrt(mysqrt);
3041 double firstterm = xc / 2 + (xc * (rcyl * rcyl -
r *
r)) / (2 * d);
3042 double secondterm = (mysqrt * yc) / (2 * d);
3043 const double x1 = firstterm + secondterm;
3044 const double x2 = firstterm - secondterm;
3045 firstterm = yc / 2 + (yc * (rcyl * rcyl -
r *
r)) / (2 * d);
3046 secondterm = (mysqrt * xc) / (2 * d);
3047 const double y1 = firstterm - secondterm;
3048 const double y2 = firstterm + secondterm;
3051 const double dist1 = (
x - x1) * (
x - x1) + (
y - y1) * (
y - y1);
3052 const double dist2 = (
x - x2) * (
x - x2) + (
y - y2) * (
y - y2);
3054 if (dist1 < dist2) {
3062 const double phi1 = std::atan2(
y - yc,
x - xc);
3065 const double delta_z =
r * deltaphi / tantheta;
3066 const double z =
z0 + delta_z;
3077 const Amg::Vector3D trackdir(cos(phidir) * sintheta, std::sin(phidir) * sintheta, costheta);
3079 const double costracksurf = std::abs(normal.unit().dot(trackdir));
3081 return std::make_pair(intersect, costracksurf);
3088 std::vector<std::pair<const Layer *, const Layer *>> &layers,
3093 std::vector<std::unique_ptr<GXFTrackState>> & states = trajectory.
trackStates();
3094 std::vector<std::unique_ptr<GXFTrackState>> oldstates = std::move(states);
3097 states.reserve(oldstates.size() + layers.size());
3105 for (
int i = 0; i <= indexoffset; i++) {
3115 for (
int i = indexoffset + 1; i < (int) oldstates.size(); i++) {
3116 const double rmeas = oldstates[i]->position().perp();
3117 const double zmeas = oldstates[i]->position().z();
3125 while (layerindex < (
int) layers.size()) {
3127 double costracksurf = 0.0;
3128 const Layer *layer =
nullptr;
3137 if (layers[layerindex].first !=
nullptr) {
3141 layer = layers[layerindex].first;
3148 if (oldstates[i]->trackParameters() !=
nullptr) {
3149 const double rlayer = cylsurf->
bounds().
r();
3150 if (std::abs(rmeas - rlayer) < std::abs(parforextrap->
position().perp() - rlayer)) {
3151 parforextrap = oldstates[i]->trackParameters();
3161 std::tie(intersect, costracksurf) =
res.value();
3167 if (cylsurf->
bounds().
r() > rmeas)
break;
3168 }
else if (layers[layerindex].second !=
nullptr) {
3174 layer = layers[layerindex].second;
3177 if (oldstates[i]->trackParameters() !=
nullptr) {
3178 const double zlayer = discsurf->
center().z();
3179 if (std::abs(zmeas - zlayer) < std::abs(parforextrap->
position().z() - zlayer)) {
3180 parforextrap = oldstates[i]->trackParameters();
3185 std::tie(intersect, costracksurf) =
res.value();
3191 if (std::abs(discsurf->
center().z()) > std::abs(zmeas))
break;
3193 throw std::logic_error(
"Unhandled surface.");
3200 const MaterialProperties *matprop = layer->layerMaterialProperties()->fullMaterial(intersect);
3201 if (matprop ==
nullptr) {
3214 const double actualx0 = X0 / costracksurf;
3215 const double de = -std::abs(
3219 (
m_p != 0.0 ? std::abs(
m_p) : std::abs(1. / currentqoverp)),
3221 const double sintheta = std::sin(parforextrap->parameters()[
Trk::theta]);
3222 const double sigmascat = std::sqrt(
m_scattool->sigmaSquare(
3224 (
m_p != 0.0 ? std::abs(
m_p) : std::abs(1. / currentqoverp)),
3228 std::unique_ptr<GXFMaterialEffects> meff = std::make_unique<GXFMaterialEffects>();
3229 meff->setDeltaE(de);
3230 meff->setScatteringSigmas(sigmascat / sintheta, sigmascat);
3231 meff->setX0(actualx0);
3232 meff->setSurface(&layer->surfaceRepresentation());
3233 meff->setMaterialProperties(matprop);
3239 std::unique_ptr<EnergyLoss> eloss;
3242 eloss = std::make_unique<EnergyLoss>(
m_elosstool->energyLoss(
3244 (
m_p != 0.0 ? std::abs(
m_p) : std::abs(1. / currentqoverp)),
3249 if (eloss !=
nullptr) {
3250 meff->setSigmaDeltaE(eloss->sigmaDeltaE());
3255 meff->setDeltaE(-5);
3257 meff->setScatteringSigmas(0, 0);
3260 meff->setSigmaDeltaE(50);
3261 if (eloss !=
nullptr) {
3262 meff->setSigmaDeltaEPos(eloss->sigmaPlusDeltaE());
3263 meff->setSigmaDeltaENeg(eloss->sigmaMinusDeltaE());
3268 "X0: " << meff->x0() <<
" qoverp: " << currentqoverp <<
3269 " sigmascat " << meff->sigmaDeltaTheta() <<
" eloss: " << meff->deltaE() <<
3270 " sigma eloss: " << meff->sigmaDeltaE()
3277 std::unique_ptr<GXFTrackState> matstate = std::make_unique<GXFTrackState>(
3279 std::unique_ptr<const TrackParameters>()
3281 matstate->setPosition(intersect);
3297 std::vector<std::pair<const Layer *, const Layer *>> & layers,
3298 std::vector<std::pair<const Layer *, const Layer *>> & upstreamlayers,
3299 const std::vector<std::unique_ptr<GXFTrackState>> & oldstates,
3308 upstreamlayers.reserve(5);
3327 const double lastz = laststate->
position().z();
3328 const double lastr = laststate->
position().perp();
3334 const double tantheta = std::tan(refpar->parameters()[
Trk::theta]);
3335 const double slope = (tantheta != 0) ? 1 / tantheta : 0;
3341 std::vector < const Layer *>::const_iterator it;
3342 std::vector < const Layer *>::const_iterator itend;
3350 it =
cache.m_posdiscs.begin();
3351 itend =
cache.m_posdiscs.end();
3353 it =
cache.m_negdiscs.begin();
3354 itend =
cache.m_negdiscs.end();
3360 for (; it != itend; ++it) {
3365 if (std::abs((*it)->surfaceRepresentation().center().z()) > std::abs(lastz)) {
3373 const DiscBounds *discbounds =
static_cast<const DiscBounds *
> (&(*it)->surfaceRepresentation().bounds());
3378 if (discbounds->
rMax() < firstr || discbounds->
rMin() > lastr) {
3382 const double rintersect = firstr + ((*it)->surfaceRepresentation().center().
z() - firstz) / slope;
3385 rintersect < discbounds->rMin() - 50 ||
3386 rintersect > discbounds->
rMax() + 50
3396 if ((*it) == endlayer) {
3408 std::abs((*it)->surfaceRepresentation().center().z()) < std::abs(firstz) ||
3411 upstreamlayers.emplace_back((
Layer *)
nullptr, (*it));
3418 (*it) != startlayer &&
3419 (std::abs((*it)->surfaceRepresentation().center().z()) > std::abs(firstz2) ||
3420 (*it) == startlayer2)
3422 layers.emplace_back((
Layer *)
nullptr, (*it));
3430 for (
const auto *barrelcylinder :
cache.m_barrelcylinders) {
3434 if (barrelcylinder->surfaceRepresentation().bounds().r() > lastr) {
3441 const double zintersect = firstz + (barrelcylinder->surfaceRepresentation().bounds().
r() - firstr) * slope;
3443 if (std::abs(zintersect - barrelcylinder->surfaceRepresentation().center().z()) >
3444 ((
const CylinderSurface*)(&barrelcylinder->surfaceRepresentation()))->bounds().halflengthZ() + 50) {
3448 if (barrelcylinder == endlayer) {
3455 if (barrelcylinder->surfaceRepresentation().bounds().r() < firstr ||
3456 barrelcylinder == startlayer) {
3457 upstreamlayers.emplace_back(barrelcylinder, (
Layer*)
nullptr);
3460 if (barrelcylinder != startlayer &&
3461 (barrelcylinder->surfaceRepresentation().bounds().r() > firstr2 ||
3462 barrelcylinder == startlayer2)) {
3463 layers.emplace_back(barrelcylinder, (
Layer*)
nullptr);
3476 const EventContext& ctx,
3487 if (!caloEntranceIsValid) {
3496 cache.m_negdiscs.empty() &&
3497 cache.m_posdiscs.empty() &&
3498 cache.m_barrelcylinders.empty()
3513 cache.m_fastmat =
false;
3530 bool hasmat =
false;
3531 int indexoffset = 0;
3532 int lastmatindex = 0;
3533 std::vector<std::unique_ptr<GXFTrackState>> & oldstates = trajectory.
trackStates();
3547 for (
int i = 0; i < (int) oldstates.size(); i++) {
3548 if (oldstates[i]->materialEffects() !=
nullptr) {
3557 if (firstsistate ==
nullptr) {
3558 if (oldstates[i]->trackParameters() ==
nullptr) {
3559 std::unique_ptr<const TrackParameters> tmppar(
m_propagator->propagateParameters(
3562 oldstates[i]->associatedSurface(),
3569 if (tmppar ==
nullptr)
return;
3571 oldstates[i]->setTrackParameters(std::move(tmppar));
3573 firstsistate = oldstates[i].get();
3575 lastsistate = oldstates[i].get();
3583 if (lastsistate ==
nullptr) {
3584 throw std::logic_error(
"No track state");
3594 std::unique_ptr<const TrackParameters> tmppar(
m_propagator->propagateParameters(
3603 if (tmppar ==
nullptr)
return;
3615 indexoffset = lastmatindex;
3630 std::vector<std::pair<const Layer *, const Layer *>> layers;
3631 std::vector<std::pair<const Layer *, const Layer *>> & upstreamlayers = trajectory.
upstreamMaterialLayers();
3645 const EventContext& ctx,
3651 if (refpar2 ==
nullptr) {
3661 std::vector<std::unique_ptr<GXFTrackState>> & states = trajectory.
trackStates();
3662 std::vector<std::unique_ptr<GXFTrackState>> matstates;
3663 std::unique_ptr< const std::vector < const TrackStateOnSurface *>,
3664 void (*)(
const std::vector<const TrackStateOnSurface *> *) >
3666 bool matvec_used=
false;
3667 std::unique_ptr<TrackParameters> startmatpar1;
3668 std::unique_ptr<TrackParameters> startmatpar2;
3679 int npseudomuon1 = 0;
3680 int npseudomuon2 = 0;
3682 for (
auto & state : states) {
3688 if (firstidhit ==
nullptr) {
3697 if (firsthit ==
nullptr) {
3698 firsthit = state->measurement();
3699 if (
cache.m_acceleration) {
3700 if (tp ==
nullptr) {
3710 if (tp ==
nullptr) {
3714 state->setTrackParameters(std::unique_ptr<const TrackParameters>(tp));
3720 lasthit = state->measurement();
3726 if (firstidhit ==
nullptr) {
3727 firstidhit = state->measurement();
3730 if ((firstidpar ==
nullptr) && (tp !=
nullptr)) {
3734 lastidhit = state->measurement();
3735 if (tp !=
nullptr) {
3740 if (firstsiliconpar ==
nullptr) {
3741 firstsiliconpar = tp;
3743 lastsiliconpar = tp;
3755 if (firstmuonhit ==
nullptr) {
3756 firstmuonhit = state->measurement();
3757 if (tp !=
nullptr) {
3761 lastmuonhit = state->measurement();
3762 if (tp !=
nullptr) {
3768 if (meff->
deltaE() == 0) {
3769 if (firstcalopar ==
nullptr) {
3770 firstcalopar = state->trackParameters();
3772 lastcalopar = state->trackParameters();
3774 if (firstmatpar ==
nullptr) {
3775 firstmatpar = state->trackParameters();
3780 std::unique_ptr<TrackParameters> refpar;
3781 AmgVector(5) newpars = refpar2->parameters();
3788 newpars[0], newpars[1], newpars[2], newpars[3], newpars[4], std::nullopt
3791 if (firstmatpar !=
nullptr) {
3796 if ((startmatpar1 ==
nullptr) || ((firstidhit !=
nullptr) && (firstmuonhit !=
nullptr))) {
3800 const double mass = trajectory.
mass();
3801 if (mass > 200 *
MeV) {
3802 const AmgVector(5) & newpars = startmatpar2->parameters();
3803 const double oldp = std::abs(1 / newpars[
Trk::qOverP]);
3806 startmatpar2 = startmatpar2->associatedSurface().createUniqueTrackParameters(
3807 newpars[0], newpars[1], newpars[2], newpars[3],
3808 sign / std::sqrt(oldp * oldp + 2 * 100 *
MeV * std::sqrt(oldp * oldp + mass * mass) + 100 *
MeV * 100 *
MeV),
3813 AmgVector(5) newpars = startmatpar1->parameters();
3816 startmatpar1 = startmatpar1->associatedSurface().createUniqueTrackParameters(
3817 newpars[0], newpars[1], newpars[2], newpars[3], newpars[4], std::nullopt
3820 newpars = startmatpar2->parameters();
3823 startmatpar2 = startmatpar2->associatedSurface().createUniqueTrackParameters(
3824 newpars[0], newpars[1], newpars[2], newpars[3], newpars[4], std::nullopt
3832 refpar->momentum().unit()
3835 const double distance = getDistance(distsol);
3838 ATH_MSG_DEBUG(
"Obtaining upstream layers from Extrapolator");
3841 std::unique_ptr<const TrackParameters> tmppar;
3843 if (firstmuonhit !=
nullptr) {
3845 if (caloEntranceIsValid) {
3848 *
cache.m_caloEntrance,
3852 if (tmppar !=
nullptr) {
3853 destsurf = &tmppar->associatedSurface();
3858 if (matvec_used)
cache.m_matTempStore.push_back( std::move(matvec) );
3863 false, matEffects) );
3866 if (matvec && !matvec->empty()) {
3867 for (
int i = (
int)matvec->size() - 1; i > -1; i--) {
3872 std::unique_ptr<GXFMaterialEffects> meff = std::make_unique<GXFMaterialEffects>(*meot);
3873 const TrackParameters * newpars = (*matvec)[i]->trackParameters() !=
nullptr ? (*matvec)[i]->trackParameters()->clone() :
nullptr;
3874 meff->setSigmaDeltaE(0);
3875 matstates.push_back(std::make_unique<GXFTrackState>(
3877 std::unique_ptr<const TrackParameters>(newpars)
3890 refpar->momentum().unit()
3893 const double distance = getDistance(distsol);
3896 ATH_MSG_DEBUG(
"Obtaining downstream ID layers from Extrapolator");
3898 std::unique_ptr<const TrackParameters> tmppar;
3899 std::unique_ptr<Surface> calosurf;
3900 if (firstmuonhit !=
nullptr) {
3902 if (caloEntranceIsValid) {
3905 *
cache.m_caloEntrance,
3910 if (tmppar !=
nullptr) {
3914 cylcalosurf =
static_cast<const CylinderSurface *
>(&tmppar->associatedSurface());
3919 disccalosurf =
static_cast<const DiscSurface *
>(&tmppar->associatedSurface());
3921 if (cylcalosurf !=
nullptr) {
3924 const double radius = cylbounds.
r();
3926 calosurf = std::make_unique<CylinderSurface>(trans, radius - 1, hlength);
3927 }
else if (disccalosurf !=
nullptr) {
3928 const double newz = (
3929 disccalosurf->
center().
z() > 0 ?
3930 disccalosurf->
center().z() - 1 :
3931 disccalosurf->
center().z() + 1
3935 disccalosurf->
center().x(),
3936 disccalosurf->
center().y(),
3941 trans.translation() << newpos;
3944 const double rmin = discbounds->
rMin();
3945 const double rmax = discbounds->
rMax();
3946 calosurf = std::make_unique<DiscSurface>(trans, rmin, rmax);
3948 destsurf = calosurf.release();
3952 if (matvec_used)
cache.m_matTempStore.push_back( std::move(matvec) );
3954 ctx, *startmatpar2, *destsurf,
alongMomentum,
false, matEffects));
3955 matvec_used =
false;
3957 if (matvec && !matvec->empty()) {
3958 for (
const auto & i : *matvec) {
3964 std::unique_ptr<GXFMaterialEffects> meff = std::make_unique<GXFMaterialEffects>(*meot);
3970 meff->setDeltaE(-5);
3973 meff->setScatteringSigmas(0, 0);
3976 meff->setSigmaDeltaE(50);
3979 const TrackParameters * newparams = i->trackParameters() !=
nullptr ? i->trackParameters()->clone() :
nullptr;
3981 matstates.push_back(std::make_unique<GXFTrackState>(
3983 std::unique_ptr<const TrackParameters>(newparams)
3995 if (
cache.m_calomat && (firstmuonhit !=
nullptr) && (firstidhit !=
nullptr)) {
3998 std::vector<MaterialEffectsOnTrack> calomeots =
m_calotool->extrapolationSurfacesAndEffects(
4007 if (calomeots.empty()) {
4012 std::unique_ptr<const TrackParameters> prevtrackpars =
unique_clone(lastidpar);
4013 if (lasthit == lastmuonhit) {
4014 for (
int i = 0; i < (int) calomeots.size(); i++) {
4017 std::unique_ptr<const TrackParameters> layerpar(
m_propagator->propagateParameters(
4020 calomeots[i].associatedSurface(),
4027 if (layerpar ==
nullptr) {
4032 std::unique_ptr<GXFMaterialEffects> meff = std::make_unique<GXFMaterialEffects>(calomeots[i]);
4035 lastcalopar = layerpar.get();
4039 const double qoverp = layerpar->parameters()[
Trk::qOverP];
4040 double qoverpbrem = 0;
4044 (firstmuonpar !=
nullptr) &&
4045 std::abs(firstmuonpar->parameters()[
Trk::qOverP]) > 1.e-9
4047 qoverpbrem = firstmuonpar->parameters()[
Trk::qOverP];
4049 const double sign = (qoverp > 0) ? 1 : -1;
4050 qoverpbrem =
sign / (1 / std::abs(qoverp) - std::abs(calomeots[i].energyLoss()->deltaE()));
4053 const AmgVector(5) & newpar = layerpar->parameters();
4055 layerpar = layerpar->associatedSurface().createUniqueTrackParameters(
4056 newpar[0], newpar[1], newpar[2], newpar[3], qoverpbrem, std::nullopt
4058 meff->setdelta_p(1000 * (qoverpbrem - qoverp));
4061 matstates.push_back(std::make_unique<GXFTrackState>(
4063 std::unique_ptr<const TrackParameters>(layerpar !=
nullptr ? layerpar->clone() :
nullptr)
4065 prevtrackpars = std::move(layerpar);
4070 firsthit == firstmuonhit &&
4071 (!
cache.m_getmaterialfromtrack || lasthit == lastidhit)
4074 for (
int i = 0; i < (int) calomeots.size(); i++) {
4076 std::unique_ptr<const TrackParameters> layerpar(
m_propagator->propagateParameters(
4079 calomeots[i].associatedSurface(),
4086 if (layerpar ==
nullptr) {
4091 std::unique_ptr<GXFMaterialEffects> meff = std::make_unique<GXFMaterialEffects>(calomeots[i]);
4100 const double qoverpbrem = layerpar->parameters()[
Trk::qOverP];
4105 (lastmuonpar !=
nullptr) &&
4106 std::abs(lastmuonpar->parameters()[
Trk::qOverP]) > 1.e-9
4110 const double sign = (qoverpbrem > 0) ? 1 : -1;
4111 qoverp =
sign / (1 / std::abs(qoverpbrem) + std::abs(calomeots[i].energyLoss()->deltaE()));
4114 meff->setdelta_p(1000 * (qoverpbrem - qoverp));
4115 const AmgVector(5) & newpar = layerpar->parameters();
4117 prevtrackpars = layerpar->associatedSurface().createUniqueTrackParameters(
4118 newpar[0], newpar[1], newpar[2], newpar[3], qoverp, std::nullopt
4122 matstates.insert(matstates.begin(), std::make_unique<GXFTrackState>(std::move(meff), std::move(layerpar)));
4127 if (lasthit == lastmuonhit &&
cache.m_extmat) {
4128 std::unique_ptr<const Trk::TrackParameters> muonpar1;
4130 if (lastcalopar !=
nullptr) {
4132 if (msEntranceIsValid) {
4136 *
cache.m_msEntrance,
4140 if (muonpar1 !=
nullptr) {
4148 rot.col(2) = trackdir;
4150 trans.linear().matrix() << rot;
4151 trans.translation() << muonpar1->position() - .1 * trackdir;
4154 std::unique_ptr<const TrackParameters> curvlinpar(
m_extrapolator->extrapolateDirectly(
4162 if (curvlinpar !=
nullptr) {
4163 muonpar1 = std::move(curvlinpar);
4167 muonpar1 = std::unique_ptr<TrackParameters>(lastcalopar->
clone());
4171 muonpar1 = std::unique_ptr<TrackParameters>(refpar->clone());
4176 if (muonpar1 !=
nullptr) {
4178 muonpar1->position(),
4179 muonpar1->momentum().unit()
4183 double distance = getDistance(distsol);
4186 0) and (firstmuonhit !=
nullptr)) {
4188 muonpar1->position(),
4189 muonpar1->momentum().unit()
4195 distance = distsol.
first();
4198 std::abs(distsol.
first()) < std::abs(distsol.
second()) ?
4204 if (distance < 0 && distsol.
numberOfSolutions() > 0 && (firstidhit ==
nullptr)) {
4205 if (firstmuonpar !=
nullptr) {
4206 AmgVector(5) newpars = firstmuonpar->parameters();
4213 newpars[0], newpars[1], newpars[2], newpars[3], newpars[4], std::nullopt
4216 std::unique_ptr<const TrackParameters> tmppar(
m_propagator->propagateParameters(
4226 if (tmppar !=
nullptr) {
4227 muonpar1 = std::move(tmppar);
4233 ATH_MSG_DEBUG(
"Obtaining downstream layers from Extrapolator");
4234 if (matvec_used)
cache.m_matTempStore.push_back( std::move(matvec) );
4237 states.back()->associatedSurface(),
4241 matvec_used =
false;
4248 if (matvec && !matvec->empty()) {
4249 for (
int j = 0; j < (int) matvec->size(); j++) {
4255 std::unique_ptr<GXFMaterialEffects> meff = std::make_unique<GXFMaterialEffects>(*meot);
4260 std::abs(meff->deltaE()) > 25 &&
4261 std::abs((*matvec)[j]->trackParameters()->position().z()) < 13000
4266 const TrackParameters * newparams = (*matvec)[j]->trackParameters() !=
nullptr ? (*matvec)[j]->trackParameters()->clone() :
nullptr;
4268 matstates.push_back(std::make_unique<GXFTrackState>(
4270 std::unique_ptr<const TrackParameters>(newparams)
4280 if (firsthit == firstmuonhit &&
cache.m_extmat && (firstcalopar !=
nullptr)) {
4281 std::unique_ptr<const Trk::TrackParameters> muonpar1;
4284 if (msEntranceIsValid) {
4288 *
cache.m_msEntrance,
4292 if (muonpar1 !=
nullptr) {
4300 rot.col(2) = trackdir;
4302 trans.linear().matrix() << rot;
4303 trans.translation() << muonpar1->position() - .1 * trackdir;
4306 std::unique_ptr<const TrackParameters> curvlinpar(
m_extrapolator->extrapolateDirectly(
4314 if (curvlinpar !=
nullptr) {
4315 muonpar1 = std::move(curvlinpar);
4319 muonpar1 = std::unique_ptr<const TrackParameters>(firstcalopar->
clone());
4325 if (muonpar1 !=
nullptr) {
4327 muonpar1->position(),
4328 muonpar1->momentum().unit()
4332 const double distance = getDistance(distsol);
4336 ATH_MSG_DEBUG(
"Collecting upstream muon material from extrapolator");
4337 if (matvec_used)
cache.m_matTempStore.push_back( std::move(matvec) );
4340 states[0]->associatedSurface(),
4344 matvec_used =
false;
4346 if (matvec && !matvec->empty()) {
4347 ATH_MSG_DEBUG(
"Retrieved " << matvec->size() <<
" material states");
4349 for (
int j = 0; j < (int) matvec->size(); j++) {
4352 if (meb !=
nullptr) {
4357 std::unique_ptr<GXFMaterialEffects> meff = std::make_unique<GXFMaterialEffects>(*meot);
4362 std::abs(meff->deltaE()) > 25 &&
4363 std::abs((*matvec)[j]->trackParameters()->position().z()) < 13000
4369 (*matvec)[j]->trackParameters() !=
nullptr
4370 ? (*matvec)[j]->trackParameters()->clone()
4373 matstates.insert(matstates.begin(), std::make_unique<GXFTrackState>(
4375 std::unique_ptr<const TrackParameters>(tmpparams)
4388 std::vector<std::unique_ptr<GXFTrackState>> & newstates = states;
4389 std::vector<std::unique_ptr<GXFTrackState>> oldstates = std::move(newstates);
4392 newstates.reserve(oldstates.size() + matstates.size());
4395 int firstlayerno = -1;
4397 if (
cache.m_acceleration) {
4401 const double cosphi = std::cos(refpar->parameters()[
Trk::phi0]);
4402 const double sinphi = std::sin(refpar->parameters()[
Trk::phi0]);
4404 for (
int i =
cache.m_acceleration ? 1 : 0; i < (
int) oldstates.size(); i++) {
4405 bool addlayer =
true;
4407 while (addlayer && layerno < (
int) matstates.size()) {
4409 const TrackParameters *layerpar = matstates[layerno]->trackParameters();
4411 const DistanceSolution distsol = oldstates[i]->associatedSurface().straightLineDistanceEstimate(
4416 const double distance = getDistance(distsol);
4424 const double trackimpact = std::abs(-refpar->position().x() * sinphi + refpar->position().y() * cosphi);
4426 if (trackimpact > cylinderradius - 5 * mm) {
4432 if (i == (
int) oldstates.size() - 1) {
4441 ATH_MSG_DEBUG(
"Material state with excessive scattering, skipping it");
4452 if (firstlayerno < 0) {
4453 firstlayerno = layerno;
4458 if ((layerpar !=
nullptr) && matEffects !=
pion && matEffects !=
muon) {
4460 const Layer *lay =
nullptr;
4462 if (tvol !=
nullptr) {
4478 if (matvec_used)
cache.m_matTempStore.push_back( std::move(matvec) );
4491 if ((persurf !=
nullptr) && (!
cache.m_acceleration || persurf->
center().perp() > 5)) {
4492 const AmgVector(5) & pars = param.parameters();
4494 pars[0], pars[1], pars[2], pars[3], pars[4], std::nullopt
4498 if (
cache.m_acceleration) {
4503 std::unique_ptr<const TrackParameters> per(
m_extrapolator->extrapolate(
4504 Gaudi::Hive::currentContext(),param, tmppersf,
oppositeMomentum,
false, matEffects));
4506 if (per ==
nullptr) {
4507 ATH_MSG_DEBUG(
"Cannot make Perigee with starting parameters");
4511 if(std::abs(per->position().z())>5000.) {
4512 ATH_MSG_WARNING(
"Pathological perigee well outside of tracking detector!! Returning nullptr");
4520 const EventContext& ctx,
4527 ATH_MSG_DEBUG(
"--> entering GlobalChi2Fitter::myfit_helper");
4542 cache.m_lastiter = 0;
4555 if (trajectory.
nDOF() < 0) {
4560 cache.m_phiweight.clear();
4561 cache.m_firstmeasurement.clear();
4562 cache.m_lastmeasurement.clear();
4565 ATH_MSG_WARNING(
"Attempt to apply material corrections with q/p=0, reject track");
4579 std::unique_ptr<const TrackParameters> per =
makePerigee(
cache, param, matEffects);
4581 if (!
cache.m_acceleration && (per ==
nullptr)) {
4591 cache.m_acceleration &&
4595 ATH_MSG_WARNING(
"Tracking Geometry Service and/or Material Updator Tool not configured");
4598 cache.m_fastmat =
false;
4603 !
cache.m_acceleration ||
4606 addMaterial(ctx,
cache, trajectory, per !=
nullptr ? per.get() : ¶m, matEffects);
4609 ctx,
cache, trajectory, per !=
nullptr ? per.get() : ¶m, matEffects);
4621 for (
const auto & state : trajectory.
trackStates()) {
4625 state->materialEffects()->deltaE() == 0
4627 scatstate2 = state.get();
4632 scatstate = state.get();
4639 if ((scatstate ==
nullptr) || (scatstate2 ==
nullptr)) {
4640 throw std::logic_error(
"Invalid scatterer");
4645 const int nstates = (int) trajectory.
trackStates().size();
4653 std::unique_ptr<const TrackParameters> nearestpar;
4654 double mindist = 99999;
4655 std::vector < GXFTrackState * >mymatvec;
4658 if ((*it).trackParameters() ==
nullptr) {
4662 const double distance = persurf
4664 (*it).trackParameters()->position(),
4665 (*it).trackParameters()->momentum().unit())
4668 const bool insideid = (
4669 (
cache.m_caloEntrance ==
nullptr) ||
4670 cache.m_caloEntrance->inside((*it).trackParameters()->position())
4674 (((*it).measurement() !=
nullptr) && insideid) || (
4675 ((*it).materialEffects() !=
nullptr) &&
4677 (*it).materialEffects()->deltaE() == 0 ||
4678 ((*it).materialEffects()->sigmaDeltaPhi() == 0 &&
4680 (*it).materialEffects()->deltaPhi() != 0
4684 const double dist = ((*it).trackParameters()->
position() -
vertex).perp();
4685 if (dist < mindist) {
4693 if (((*it).materialEffects() !=
nullptr) && distance > 0) {
4694 mymatvec.push_back(it.get());
4698 if (nearestpar ==
nullptr) {
4702 for (
auto & state : mymatvec) {
4704 const Surface &matsurf = state->associatedSurface();
4706 nearestpar->position(), nearestpar->momentum().unit());
4708 const double distance = getDistance(distsol);
4714 std::unique_ptr<const TrackParameters> tmppar(
m_propagator->propagateParameters(
4724 if (tmppar ==
nullptr) {
4736 if (tmppar ==
nullptr) {
4746 AmgVector(5) newpars = tmppar->parameters();
4748 if (state->materialEffects()->sigmaDeltaE() > 0) {
4749 newpars[
Trk::qOverP] += .001 * state->materialEffects()->delta_p();
4752 const double de = std::abs(state->materialEffects()->deltaE());
4753 const double oldp = std::abs(1 / newpars[
Trk::qOverP]);
4754 const double newp2 = oldp * oldp - 2 * de * std::sqrt(mass * mass + oldp * oldp) + de * de;
4760 nearestpar = tmppar->associatedSurface().createUniqueTrackParameters(
4761 newpars[0], newpars[1], newpars[2], newpars[3], newpars[4], std::nullopt
4765 std::unique_ptr<Trk::TrackParameters> tmpPars(
m_propagator->propagateParameters(
4776 if (tmpPars !=
nullptr) {
4777 per.reset(
static_cast < const Perigee *
>(tmpPars.release()));
4782 const double oldp = 1. / std::abs(per->parameters()[
Trk::qOverP]);
4784 const double newp = std::sqrt(oldp * oldp + 2 * toteloss * std::sqrt(oldp * oldp + mass * mass) + toteloss * toteloss);
4785 AmgVector(5) params = per->parameters();
4788 per = per->associatedSurface().createUniqueTrackParameters(
4789 params[0], params[1], params[2], params[3], params[4], std::nullopt
4793 if (per ==
nullptr) {
4810 }
else if (per ==
nullptr) {
4827 const AmgVector(5) & pars = per->parameters();
4828 per = per->associatedSurface().createUniqueTrackParameters(
4829 pars[0], pars[1], pars[2], pars[3], 0, std::nullopt
4835 if (per !=
nullptr) {
4845 Eigen::MatrixXd a_inv;
4846 a.resize(nfitpar, nfitpar);
4851 derivPool.setZero();
4853 for (std::unique_ptr<GXFTrackState> & state : trajectory.
trackStates()) {
4854 if (state->materialEffects() !=
nullptr) {
4857 state->setDerivatives(derivPool);
4860 bool doderiv =
true;
4861 const int tmpminiter =
cache.m_miniter;
4863 for (
int it = 0; it <
m_maxit; ++it) {
4864 cache.m_lastiter = it;
4870 cache.m_miniter = tmpminiter;
4875 cache.m_fittercode =
4885 cache.m_miniter = tmpminiter;
4892 const double redchi2 = (trajectory.
nDOF() > 0) ? trajectory.
chi2() / trajectory.
nDOF() : 0;
4893 const double prevredchi2 = (trajectory.
nDOF() > 0) ? trajectory.
prevchi2() / trajectory.
nDOF() : 0;
4902 (redchi2 < prevredchi2 &&
4903 (redchi2 > prevredchi2 - 1 || redchi2 < 2)) ||
4904 nsihits + ntrthits == nhits
4909 if (it != 1 || nsihits != 0 || trajectory.
nDOF() <= 0 || trajectory.
chi2() / trajectory.
nDOF() <= 3) {
4914 cache.m_miniter = tmpminiter;
4924 if (ntrtprechits+ntrttubehits) {
4925 phf = float(ntrtprechits)/float(ntrtprechits+ntrttubehits);
4927 if (phf<m_minphfcut && it>=3) {
4928 if ((ntrtprechits+ntrttubehits)>=15) {
4932 ATH_MSG_DEBUG(
"Iter = " << it <<
" | nTRTStates = " << ntrthits
4933 <<
" | nTRTPrecHits = " << ntrtprechits
4934 <<
" | nTRTTubeHits = " << ntrttubehits
4935 <<
" | nOutliers = "
4944 cache.m_miniter = tmpminiter;
4953 cache.m_miniter = tmpminiter;
4955 if (trajectory.
prefit() == 0) {
4959 Eigen::LLT < Eigen::MatrixXd >
const lltOfW(
a);
4960 if (lltOfW.info() == Eigen::Success) {
4964 const int ncols =
a.cols();
4965 Amg::MatrixX const weightInvAMG = Amg::MatrixX::Identity(ncols, ncols);
4966 a_inv = lltOfW.solve(weightInvAMG);
4977 (runOutlier ||
cache.m_sirecal) &&
4985 if (traj != &trajectory) {
4990 finaltrajectory = traj;
4997 if (!
cache.m_acceleration && (finaltrajectory->
prefit() == 0)) {
4998 if (nperpars == 5) {
4999 for (
int i = 0; i <
a.cols(); i++) {
5000 a_inv(4, i) *= .001;
5001 a_inv(i, 4) *= .001;
5005 int scatterPos = nperpars + 2 * nscat;
5006 for (
int bremno = 0; bremno < nbrem; bremno++, scatterPos++) {
5007 for (
int i = 0; i <
a.cols(); i++) {
5008 a_inv(scatterPos, i) *= .001;
5009 a_inv(i, scatterPos) *= .001;
5016 for (
int i = 0; i < nperparams; i++) {
5017 for (
int j = 0; j < nperparams; j++) {
5018 (errmat) (j, i) = a_inv(j, i);
5023 (errmat) (4, 4) = 1e-20;
5027 std::unique_ptr<const TrackParameters> measper(
5029 perpars[0], perpars[1], perpars[2], perpars[3], perpars[4], std::move(errmat)
5035 cache.m_fullcovmat = a_inv;
5039 std::unique_ptr<Track> track =
nullptr;
5041 if (finaltrajectory->
prefit() > 0) {
5042 if (finaltrajectory != &trajectory) {
5044 delete finaltrajectory;
5063 (track !=
nullptr) && (
5064 track->fitQuality()->numberDoF() != 0 &&
5065 track->fitQuality()->chiSquared() / track->fitQuality()->numberDoF() > cut
5068 track.reset(
nullptr);
5072 if (track ==
nullptr) {
5076 if (finaltrajectory != &trajectory) {
5077 delete finaltrajectory;
5080 return track.release();
5084 const EventContext& ctx,
5089 int & bremno_maxbrempull,
5094 std::vector<std::unique_ptr<GXFTrackState>> & states = trajectory.
trackStates();
5114 const int nmeas = (int)
res.size();
5124 const int nDOF = trajectory.
nDOF();
5125 const bool doNewPseudoMeasurements = (
5129 std::abs((trajectory.
prevchi2() - trajectory.
chi2()) / nDOF) < 15 &&
5131 (nperpars == 0 || nidhits > 0)
5142 double maxbrempull = -0.2;
5151 for (
int hitno = 0; hitno < (int) states.size(); hitno++) {
5152 std::unique_ptr<GXFTrackState> & state = states[hitno];
5166 doNewPseudoMeasurements &&
5168 !state->associatedSurface().isFree() &&
5169 !state->isRecalibrated()
5172 covMatrix(0, 0) = 100;
5174 std::unique_ptr<const PseudoMeasurementOnTrack> newpseudo = std::make_unique<const PseudoMeasurementOnTrack>(
5176 std::move(covMatrix),
5180 state->setMeasurement(std::move(newpseudo));
5181 measbase = state->measurement();
5188 double *errors = state->measurementErrors();
5190 for (
int i = 0; i < 5; i++) {
5210 res[measno] = residuals[i];
5216 res[measno] = -std::remainder(-
res[measno], 2 *
M_PI);
5225 double *errors = state->measurementErrors();
5226 for (
int i = 0; i < 5; i++) {
5227 if (errors[i] > 0) {
5228 error[measno] = errors[i];
5239 ((trajectory.
prefit() == 0) || state->materialEffects()->deltaE() == 0)
5241 const double deltaPhi = state->materialEffects()->deltaPhi();
5242 const double measDeltaPhi = state->materialEffects()->measuredDeltaPhi();
5243 const double sigma2deltaPhi = std::pow(state->materialEffects()->sigmaDeltaPhi(), 2);
5244 const double deltaTheta = state->materialEffects()->deltaTheta();
5245 const double sigma2deltaTheta = std::pow(state->materialEffects()->sigmaDeltaTheta(), 2);
5247 if (trajectory.
prefit() != 1) {
5248 b[nperpars + 2 * scatno] -= (
deltaPhi - measDeltaPhi) / sigma2deltaPhi;
5249 b[nperpars + 2 * scatno + 1] -= deltaTheta / sigma2deltaTheta;
5251 b[nperpars + scatno] -= deltaTheta / sigma2deltaTheta;
5256 deltaTheta * deltaTheta / sigma2deltaTheta
5265 if ((state->materialEffects() !=
nullptr) && state->materialEffects()->sigmaDeltaE() > 0) {
5266 double averagenergyloss = std::abs(state->materialEffects()->deltaE());
5267 const double qoverpbrem = limitInversePValue(1000 * states[hitno]->trackParameters()->parameters()[
Trk::qOverP]);
5268 const double qoverp = limitInversePValue(qoverpbrem - state->materialEffects()->delta_p());
5269 const double pbrem = 1. / std::abs(qoverpbrem);
5270 const double p = 1. / std::abs(qoverp);
5271 const double mass = .001 * trajectory.
mass();
5272 const double energy = std::sqrt(p * p + mass * mass);
5273 const double bremEnergy = std::sqrt(pbrem * pbrem + mass * mass);
5275 const double resMaterial = .001 * averagenergyloss - energy + bremEnergy;
5276 res[nmeas - nbrem + bremno] = resMaterial;
5278 const double sigde = state->materialEffects()->sigmaDeltaE();
5279 const double sigdepos = state->materialEffects()->sigmaDeltaEPos();
5280 const double sigdeneg = state->materialEffects()->sigmaDeltaENeg();
5282 double errorMaterial = .001 * state->materialEffects()->sigmaDeltaE();
5283 error[nmeas - nbrem + bremno] = errorMaterial;
5294 if (state->materialEffects()->isKink()) {
5295 maxbrempull = -999999999;
5296 state_maxbrempull =
nullptr;
5300 cache.m_asymeloss &&
5302 trajectory.
prefit() == 0 &&
5304 sigde != sigdepos &&
5307 const double elosspull = resMaterial / errorMaterial;
5309 if (trajectory.
mass() > 100) {
5314 if (std::abs(elosspull) > 1) {
5315 if (elosspull < -1) {
5316 state->materialEffects()->setSigmaDeltaE(sigdepos);
5318 state->materialEffects()->setSigmaDeltaE(sigdeneg);
5321 errorMaterial = .001 * state->materialEffects()->sigmaDeltaE();
5322 error[nmeas - nbrem + bremno] = errorMaterial;
5334 !state->materialEffects()->isKink() && (
5335 (
m_fixbrem == -1 && elosspull < maxbrempull) ||
5339 bremno_maxbrempull = bremno;
5340 state_maxbrempull = state.get();
5341 maxbrempull = elosspull;
5350 trajectory.
prefit() == 0 &&
5351 state->materialEffects()->sigmaDeltaPhi() == 0 &&
5352 state->materialEffects()->isMeasuredEloss() &&
5353 resMaterial / (.001 * state->materialEffects()->sigmaDeltaEAve()) > 2.5
5355 const TrackParameters* parforcalo = states[hitno - 2]->trackParameters();
5358 std::vector<MaterialEffectsOnTrack> calomeots =
5371 if (calomeots.size() == 3) {
5372 averagenergyloss = std::abs(calomeots[1].energyLoss()->deltaE());
5373 const double newres = .001 * averagenergyloss - energy + bremEnergy;
5374 const double newerr = .001 * calomeots[1].energyLoss()->sigmaDeltaE();
5376 const double oldPull = resMaterial / errorMaterial;
5377 const double newPull = newres / newerr;
5379 if (std::abs(newPull) < std::abs(oldPull)) {
5380 ATH_MSG_DEBUG(
"Changing from measured to parametrized energy loss");
5382 state->materialEffects()->setEloss(std::unique_ptr<EnergyLoss>(calomeots[1].energyLoss()->clone()));
5383 state->materialEffects()->setSigmaDeltaE(calomeots[1].energyLoss()->sigmaDeltaE());
5384 res[nmeas - nbrem + bremno] = newres;
5385 error[nmeas - nbrem + bremno] = newerr;
5389 state->materialEffects()->setMeasuredEloss(
false);
5399 for (
int imeas = 0; imeas < nmeas; imeas++) {
5400 if (
error[imeas] == 0) {
5421 const double oldChi2 = trajectory.
prevchi2();
5422 const double newChi2 = trajectory.
chi2();
5427 const double nDOF = trajectory.
nDOF();
5428 const double oldRedChi2 = (nDOF > 0) ? oldChi2 / nDOF : 0;
5429 const double newRedChi2 = (nDOF > 0) ? newChi2 / nDOF : 0;
5432 trajectory.
prefit() > 0 && (
5433 (newRedChi2 < 2 && it != 0) ||
5434 (newRedChi2 < oldRedChi2 + .1 && std::abs(newRedChi2 - oldRedChi2) < 1 && it != 1)
5447 int miniter = (nsihits != 0 && nsihits + ntrthits == nhits) ? 1 : 2;
5448 miniter = std::max(miniter,
cache.m_miniter);
5450 if (it >= miniter && std::abs(oldChi2 - newChi2) < 1) {
5457 const int bremno_maxbrempull,
5463 if (state_maxbrempull ==
nullptr) {
5475 const int nmeas = (int)
res.size();
5478 const double oldError =
error[nmeas - nbrem + bremno_maxbrempull];
5480 error[nmeas - nbrem + bremno_maxbrempull] = newError;
5483 if (
a.cols() != nFitPars) {
5487 const double errorRatio = oldError / newError;
5488 const double errorReductionRatio = 1 - std::pow(errorRatio, 2);
5491 for (
int i = 0; i < nFitPars; i++) {
5492 if (weightderiv(nmeas - nbrem + bremno_maxbrempull, i) == 0) {
5496 for (
int j = i; j < nFitPars; j++) {
5497 const double newaij =
a(i, j) - errorReductionRatio *
5498 weightderiv(nmeas - nbrem + bremno_maxbrempull, i) *
5499 weightderiv(nmeas - nbrem + bremno_maxbrempull, j);
5501 a.fillSymmetric(i, j, newaij);
5503 weightderiv(nmeas - nbrem + bremno_maxbrempull, i) *= errorRatio;
5512 std::vector<std::unique_ptr<GXFTrackState>> & states = trajectory.
trackStates();
5525 const int nmeas = (int) weightderiv.rows();
5527 for (std::unique_ptr<GXFTrackState> & state : states) {
5531 const auto [scatmin, scatmax] = std::minmax(scatno, nscatupstream);
5532 const auto [bremmin, bremmax] = std::minmax(bremno, nbremupstream);
5539 const double sinStereo =
5541 state->sinStereo() :
5543 const double cosStereo =
5545 std::sqrt(1 - std::pow(sinStereo, 2)) :
5553 auto getThisDeriv = [sinStereo, cosStereo, &derivatives](
int i,
int j) ->
double {
5554 if (i == 0 && sinStereo != 0) {
5555 return derivatives(0, j) * cosStereo + sinStereo * derivatives(1, j);
5557 return derivatives(i, j);
5561 for (
int i = 0; i < 5; i++) {
5576 if (i == 0 && sinStereo != 0) {
5577 weightderiv.row(measno).head(cols) =
5578 (derivatives.row(0).
head(cols) * cosStereo +
5579 sinStereo * derivatives.row(1).head(cols)) /
5582 weightderiv.row(measno).head(cols) = derivatives.row(i).head(cols) /
error[measno];
5586 for (
int j = scatmin; j < scatmax; j++) {
5587 if (trajectory.
prefit() == 1) {
5588 const int index = nperparams + j;
5591 const int index = nperparams + 2 * j;
5593 weightderiv(measno,
index + 1) = getThisDeriv(i,
index + 1) /
error[measno];
5597 for (
int j = bremmin; j < bremmax; j++) {
5598 const int index = j + nperparams + 2 * nscat;
5605 double *errors = state->measurementErrors();
5606 for (
int i = 0; i < 5; i++) {
5607 if (errors[i] > 0) {
5613 ((trajectory.
prefit() == 0) || state->materialEffects()->deltaE() == 0)
5618 if ((state->materialEffects() !=
nullptr) && state->materialEffects()->sigmaDeltaE() > 0) {
5620 const double qoverpbrem = limitInversePValue(1000 * state->trackParameters()->parameters()[
Trk::qOverP]);
5621 const double qoverp = limitInversePValue(qoverpbrem - state->materialEffects()->delta_p());
5623 const double mass = .001 * trajectory.
mass();
5625 const auto thisMeasurementIdx{nmeas - nbrem + bremno};
5629 auto multiplier = [] (
double mass,
double qOverP){
5632 const auto qoverpTerm {multiplier(mass, qoverp) /
error[thisMeasurementIdx]};
5633 const auto qoverpBremTerm {multiplier(mass, qoverpbrem) /
error[thisMeasurementIdx]};
5636 weightderiv(thisMeasurementIdx, 4) = qoverpBremTerm - qoverpTerm;
5639 const auto bremNoBase = nperparams + 2 * nscat;
5640 if (bremno < nbremupstream) {
5641 weightderiv(thisMeasurementIdx, bremNoBase + bremno) = qoverpTerm;
5642 for (
int bremno2 = bremno + 1; bremno2 < nbremupstream; bremno2++) {
5643 weightderiv(thisMeasurementIdx, bremNoBase + bremno2) = qoverpTerm - qoverpBremTerm;
5646 weightderiv(thisMeasurementIdx, bremNoBase + bremno) = qoverpBremTerm;
5647 for (
int bremno2 = nbremupstream; bremno2 < bremno; bremno2++) {
5648 weightderiv(thisMeasurementIdx, bremNoBase + bremno2) = qoverpBremTerm - qoverpTerm;
5667 const int nMeas = (int)
res.size();
5669 cache.m_firstmeasurement.resize(nFitPars);
5670 cache.m_lastmeasurement.resize(nFitPars);
5672 for (
int i = 0; i < nPerPars; i++) {
5673 cache.m_firstmeasurement[i] = 0;
5674 cache.m_lastmeasurement[i] = nMeas - nBrem;
5681 const std::unique_ptr<GXFTrackState> & state = trajectory.
trackStates()[i];
5684 if (meff ==
nullptr) {
5685 measno += state->numberOfMeasuredParameters();
5689 const int firstMeasurement = i < nUpstreamStates ? 0 : measno;
5690 const int lastMeasurement = i < nUpstreamStates ? measno : nMeas - nBrem;
5693 && (trajectory.
prefit() == 0 || meff->
deltaE() == 0)) {
5694 const int scatterPos = nPerPars + 2 * scatno;
5696 cache.m_firstmeasurement[scatterPos] = firstMeasurement;
5697 cache.m_lastmeasurement[scatterPos] = lastMeasurement;
5699 cache.m_firstmeasurement[scatterPos + 1] = firstMeasurement;
5700 cache.m_lastmeasurement[scatterPos + 1] = lastMeasurement;
5706 const int bremPos = nPerPars + nScatPars + bremno;
5708 cache.m_firstmeasurement[bremPos] = firstMeasurement;
5709 cache.m_lastmeasurement[bremPos] = lastMeasurement;
5730 const int nMeas = (int)
res.size();
5732 for (
int k = 0; k < nFitPars; k++) {
5733 const int minMeasK =
cache.m_firstmeasurement[k];
5734 const int maxMeasK =
cache.m_lastmeasurement[k];
5741 for (
int measno = minMeasK; measno < maxMeasK; measno++) {
5742 b[k] +=
res[measno] * (1. /
error[measno]) * weightDeriv(measno, k);
5752 if (k == 4 || k >= nPerPars + nScatPars) {
5753 for (
int measno = nMeas - nBrem; measno < nMeas; measno++) {
5754 b[k] +=
res[measno] * (1. /
error[measno]) * weightDeriv(measno, k);
5768 for (
int k = 0; k < nFitPars; k++) {
5769 for (
int l = k; l < nFitPars; l++) {
5770 const int minMeas = std::max(
cache.m_firstmeasurement[k],
cache.m_firstmeasurement[l]);
5771 const int maxMeas = std::min(
cache.m_lastmeasurement[k],
cache.m_lastmeasurement[l]);
5774 for (
int measno = minMeas; measno < maxMeas; measno++) {
5775 a_kl += weightDeriv(measno, k) * weightDeriv(measno, l);
5778 a.fillSymmetric(l, k, a_kl);
5796 const int nMeas = (int)
res.size();
5803 for (
int k = nPerPars; k < nPerPars + nScatPars; k += 2) {
5804 a(k, k) += 1. / std::pow(scatSigmas[scatno].first, 2);
5805 a(k + 1, k + 1) += 1. / std::pow(scatSigmas[scatno].second, 2);
5813 for (
int measno = nMeas - nBrem; measno < nMeas; measno++) {
5814 for (
int k = 4; k < nFitPars; k++) {
5816 k = nPerPars + nScatPars;
5819 for (
int l = k; l < nFitPars; l++) {
5821 l = nPerPars + nScatPars;
5824 const double a_kl =
a(l, k) + weightDeriv(measno, k) * weightDeriv(measno, l);
5825 a.fillSymmetric(l, k, a_kl);
5837 const double oldRedChi2,
5838 const double newRedChi2
5846 bool weightChanged =
false;
5853 double newPhiWeight = 1.1;
5854 double newThetaWeight = 1.001;
5855 if (trajectory.
prefit() == 0) {
5861 newPhiWeight = 1.00000001;
5862 }
else if (it == 1) {
5863 newPhiWeight = 1.0000001;
5864 }
else if (it <= 3) {
5865 newPhiWeight = 1.0001;
5866 }
else if (it <= 6) {
5867 newPhiWeight = 1.01;
5870 if (newRedChi2 > oldRedChi2 - 1 && newRedChi2 < oldRedChi2) {
5871 newPhiWeight = 1.0001;
5872 newThetaWeight = 1.0001;
5873 }
else if (newRedChi2 > oldRedChi2 - 25 && newRedChi2 < oldRedChi2) {
5874 newPhiWeight = 1.001;
5875 newThetaWeight = 1.0001;
5882 std::size_t scatno = 0;
5887 for (
const auto & state : trajectory.
trackStates()) {
5890 if (meff ==
nullptr) {
5894 const bool isValidPlaneSurface =
5896 static_cast<const PlaneSurface *
>(&state->associatedSurface()) !=
nullptr;
5901 if (meff->
deltaE() == 0 || (trajectory.
prefit() == 0 && isValidPlaneSurface)) {
5902 weightChanged =
true;
5904 const int scatNoIndex = 2 * scatno + nPerPars;
5907 if (scatno >=
cache.m_phiweight.size()) {
5908 std::stringstream message;
5909 message <<
"scatno is out of range " << scatno <<
" !< " <<
cache.m_phiweight.size();
5910 throw std::range_error(message.str());
5918 a(scatNoIndex, scatNoIndex) /=
cache.m_phiweight[scatno];
5921 cache.m_phiweight[scatno] = newPhiWeight;
5922 a(scatNoIndex, scatNoIndex) *= newPhiWeight;
5923 }
else if (trajectory.
prefit() >= 2) {
5924 a(scatNoIndex, scatNoIndex) *= newPhiWeight;
5925 a(scatNoIndex + 1, scatNoIndex + 1) *= newThetaWeight;
5952 trajectory.
prefit() == 2 &&
5955 (newRedChi2 < oldRedChi2 - 25 || newRedChi2 > oldRedChi2)
5960 return weightChanged;
5969 std::size_t scatno = 0;
5978 if (scatno >=
cache.m_phiweight.size()) {
5979 std::stringstream message;
5980 message <<
"scatno is out of range " << scatno <<
" !< " <<
cache.m_phiweight.size();
5981 throw std::range_error(message.str());
5984 const bool isValidPlaneSurface =
5986 static_cast<const PlaneSurface *
>(&state->associatedSurface()) !=
nullptr;
5988 if (meff->
deltaE() == 0 || isValidPlaneSurface) {
5989 const int scatNoIndex = 2 * scatno + nPerPars;
5990 a(scatNoIndex, scatNoIndex) /=
cache.m_phiweight[scatno];
5991 cache.m_phiweight[scatno] = 1;
6005 const EventContext& ctx,
6014 const int nDOFold = trajectory.
nDOF();
6015 const double oldChi2 = trajectory.
chi2();
6016 const double oldRedChi2 = nDOFold > 0 ? oldChi2 / nDOFold : 0;
6018 if (
cache.m_phiweight.empty()) {
6038 int bremno_maxbrempull = 0;
6054 if ((state_maxbrempull !=
nullptr) && trajectory.
converged()) {
6063 const int nDOFnew = trajectory.
nDOF();
6064 const double newChi2 = trajectory.
chi2();
6065 const double newRedChi2 = nDOFnew > 0 ? newChi2 / nDOFnew : 0;
6067 ATH_MSG_DEBUG(
"old chi2: " << oldChi2 <<
"/" << nDOFold <<
"=" << oldRedChi2 <<
6068 ", new chi2: " << newChi2 <<
"/" << nDOFnew <<
"=" << newRedChi2);
6079 if (
cache.m_firstmeasurement.empty()) {
6103 if (doDeriv || weightChanged) {
6114 if (trajectory.
prefit() == 0) {
6120 if (nSiHits + nTrtHits !=
nHits) {
6127 (newRedChi2 < 2 || (newRedChi2 < oldRedChi2 && newRedChi2 > oldRedChi2 - .5))
6149 Eigen::LLT<Eigen::MatrixXd>
const llt(lu_m);
6151 if (llt.info() != Eigen::Success) {
6175 double d0 = refpar->parameters()[
Trk::d0];
6176 double z0 = refpar->parameters()[
Trk::z0];
6179 double qoverp = refpar->parameters()[
Trk::qOverP];
6181 if (nperparams > 0) {
6182 d0 += deltaParameters[0];
6183 z0 += deltaParameters[1];
6184 phi += deltaParameters[2];
6185 theta += deltaParameters[3];
6186 qoverp = (trajectory.
m_straightline) ? 0 : .001 * deltaParameters[4] + qoverp;
6198 std::vector < std::pair < double, double >>&scatangles = trajectory.
scatteringAngles();
6199 for (
int i = 0; i < nscat; i++) {
6200 scatangles[i].first += deltaParameters[2 * i + nperparams];
6201 scatangles[i].second += deltaParameters[2 * i + nperparams + 1];
6207 std::vector < double >&delta_ps = trajectory.
brems();
6208 for (
int i = 0; i < nbrem; i++) {
6209 delta_ps[i] += deltaParameters[nperparams + 2 * nscat + i];
6215 std::unique_ptr<const TrackParameters> newper(
6235 const EventContext& evtctx
6246 if (!splitProbContainer.
isValid()) {
6250 std::vector<std::unique_ptr<GXFTrackState>> & states = trajectory.
trackStates();
6257 for (
size_t stateno = 0; stateno < states.size(); stateno++) {
6262 measno += states[stateno-1]->numberOfMeasuredParameters();
6265 std::unique_ptr<GXFTrackState> & state = states[stateno];
6277 const auto *
const rot =
static_cast<const RIO_OnTrack *
>(pMeas);
6278 prd = rot->prepRawData();
6288 const auto &splitProb = splitProbContainer->splitProbability(
pixelCluster);
6289 if (!splitProb.isSplit()) {
6290 ATH_MSG_DEBUG(
"Pixel cluster is not split so no need to update" );
6294 std::unique_ptr < const RIO_OnTrack > newrot;
6295 double *olderror = state->measurementErrors();
6298 double newerror[5] = {-1,-1,-1,-1,-1};
6299 double newres[2] = {-1,-1};
6301 newrot.reset(
m_ROTcreator->correct(*prd, *trackpars, evtctx));
6306 const Amg::MatrixX & covmat = newrot->localCovariance();
6308 newerror[0] = std::sqrt(covmat(0, 0));
6309 newres[0] = newrot->localParameters()[
Trk::locX] - trackpars->parameters()[
Trk::locX];
6310 newerror[1] = std::sqrt(covmat(1, 1));
6311 newres[1] = newrot->localParameters()[
Trk::locY] - trackpars->parameters()[
Trk::locY];
6313 if (
a.cols() != nfitpars) {
6318 for(
int k =0; k<2; k++ ){
6319 const double oldres =
res[measno+k];
6320 res[measno+k] = newres[k];
6321 err[measno+k] = newerror[k];
6323 for (
int i = 0; i < nfitpars; i++) {
6324 if (weightderiv(measno+k, i) == 0) {
6328 b[i] -= weightderiv(measno+k, i) * (oldres / olderror[k] - (newres[k] * olderror[k]) / (newerror[k] * newerror[k]));
6330 for (
int j = i; j < nfitpars; j++) {
6334 weightderiv(measno+k, i) *
6335 weightderiv(measno+k, j) *
6336 ((olderror[k] * olderror[k]) / (newerror[k] * newerror[k]) - 1)
6340 weightderiv(measno+k, i) *= olderror[k] / newerror[k];
6344 state->setMeasurement(std::move(newrot));
6345 state->setMeasurementErrors(newerror);
6360 const EventContext& ctx
6368 std::vector<std::unique_ptr<GXFTrackState>> & states = trajectory.
trackStates();
6374 if (
a.cols() != nfitpars) {
6382 bool outlierremoved =
false;
6383 bool hitrecalibrated =
false;
6385 for (
int stateno = 0; stateno < (int) states.size(); stateno++) {
6386 std::unique_ptr<GXFTrackState> & state = states[stateno];
6394 std::abs(state->trackParameters()->parameters()[
Trk::driftRadius]) > 1.05 * state->associatedSurface().bounds().r()
6399 outlierremoved =
true;
6401 double *errors = state->measurementErrors();
6402 const double olderror = errors[0];
6406 for (
int i = 0; i < nfitpars; i++) {
6407 if (weightderiv(measno, i) == 0) {
6411 b[i] -=
res[measno] * weightderiv(measno, i) / olderror;
6413 for (
int j = i; j < nfitpars; j++) {
6416 a(i, j) - weightderiv(measno, i) * weightderiv(measno, j)
6419 weightderiv(measno, i) = 0;
6423 }
else if (trtrecal) {
6424 double *errors = state->measurementErrors();
6425 const double olderror = errors[0];
6427 const auto *
const thisMeasurement{state->measurement()};
6436 const double trackradius = state->trackParameters()->parameters()[
Trk::driftRadius];
6438 std::unique_ptr<const Trk::RIO_OnTrack> newrot =
nullptr;
6439 const double distance = std::abs(std::abs(trackradius) - dcradius);
6441 if (distance < scalefactor * dcerror && (olderror > 1. || trackradius * oldradius < 0)) {
6443 }
else if (distance > scalefactor * dcerror && olderror < 1.) {
6447 if (newrot !=
nullptr) {
6449 hitrecalibrated =
true;
6453 if ((measno < 0) or (measno >= (
int)
res.size())) {
6454 throw std::runtime_error(
6455 "'res' array index out of range in TrkGlobalChi2Fitter/src/GlobalChi2Fitter.cxx:" + std::to_string(__LINE__)
6459 const double oldres =
res[measno];
6460 const double newres = newradius - state->trackParameters()->parameters()[
Trk::driftRadius];
6461 errors[0] = newerror;
6462 state->setMeasurement(std::move(newrot));
6466 for (
int i = 0; i < nfitpars; i++) {
6467 if (weightderiv(measno, i) == 0) {
6471 b[i] -= weightderiv(measno, i) * (oldres / olderror - (newres * olderror) / (newerror * newerror));
6473 for (
int j = i; j < nfitpars; j++) {
6477 !
cache.m_phiweight.empty() &&
6480 i < nperpars + 2 * nscats &&
6481 (i - nperpars) % 2 == 0
6483 weight =
cache.m_phiweight[(i - nperpars) / 2];
6488 a(i, j) + weightderiv(measno, i) * weightderiv(measno, j) * ((olderror * olderror) / (newerror * newerror) - 1) * weight
6491 weightderiv(measno, i) *= olderror / newerror;
6494 res[measno] = newres;
6495 err[measno] = newerror;
6504 measno += state->numberOfMeasuredParameters();
6508 if (trajectory.
nDOF() < 0) {
6513 if (outlierremoved || hitrecalibrated) {
6517 cache.m_miniter = it + 2;
6522 const EventContext& ctx,
6530 bool trackok =
false;
6532 std::unique_ptr < GXFTrajectory > cleanup_oldtrajectory;
6534 std::unique_ptr < GXFTrajectory > cleanup_newtrajectory;
6540 while (!trackok && oldtrajectory->
nDOF() > 0) {
6542 std::vector<std::unique_ptr<GXFTrackState>> & states = oldtrajectory->
trackStates();
6550 if (nhits != nsihits) {
6554 double maxsipull = -1;
6556 int hitno_maxsipull = -1;
6557 int measno_maxsipull = -1;
6558 int stateno_maxsipull = 0;
6570 for (
int stateno = 0; stateno < (int) states.size(); stateno++) {
6571 std::unique_ptr<GXFTrackState> & state = states[stateno];
6577 double *errors = state->measurementErrors();
6579 const Amg::MatrixX & hitcov = state->measurement()->localCovariance();
6580 const double sinstereo = state->sinStereo();
6581 const double cosstereo = (sinstereo == 0) ? 1 : std::sqrt(1 - sinstereo * sinstereo);
6582 double weight1 = -1;
6584 if (hitcov(0, 0) > trackcov(0, 0)) {
6585 if (sinstereo == 0) {
6586 weight1 = errors[0] * errors[0] - trackcov(0, 0);
6588 weight1 = errors[0] * errors[0] - (
6589 trackcov(0, 0) * cosstereo * cosstereo + 2 *
6590 trackcov(1, 0) * cosstereo * sinstereo + trackcov(1, 1) * sinstereo * sinstereo
6595 const double weight2 = (
6597 errors[1] * errors[1] - trackcov(1, 1) :
6601 double sipull1 = weight1 > 0 ? std::abs(
res[measno] / std::sqrt(weight1)) : -1;
6602 const double sipull2 = (
6604 std::abs(
res[measno + 1] / std::sqrt(weight2)) :
6607 sipull1 = std::max(sipull1, sipull2);
6609 if (sipull1 > maxsipull) {
6610 maxsipull = sipull1;
6611 measno_maxsipull = measno;
6612 state_maxsipull = state.get();
6613 stateno_maxsipull = stateno;
6614 hitno_maxsipull = hitno;
6625 measno += state->numberOfMeasuredParameters();
6629 const double maxpull = maxsipull;
6631 ATH_MSG_DEBUG(
" maxsipull: " << maxsipull <<
" hitno_maxsipull: " <<
6632 hitno_maxsipull <<
" n3sigma: " << n3sigma <<
" cut: " << cut <<
" cut2: " << cut2);
6643 state_maxsipull = oldtrajectory->
trackStates()[stateno_maxsipull].get();
6646 const auto *
const rot =
static_cast<const RIO_OnTrack *
>(pMeas);
6647 prd = rot->prepRawData();
6649 std::unique_ptr < const RIO_OnTrack > broadrot;
6654 Amg::VectorX parameterVector = trackpar_maxsipull->parameters();
6655 const std::unique_ptr<const TrackParameters> trackparForCorrect(
6664 state_maxsipull->trackCovariance())
6668 newerror[0] = newerror[1] = newerror[2] = newerror[3] = newerror[4] = -1;
6669 double newpull = -1;
6670 double newpull1 = -1;
6671 double newpull2 = -1;
6672 double newres1 = -1;
6673 double newres2 = -1;
6674 double newsinstereo = 0;
6687 const Amg::MatrixX & covmat = broadrot->localCovariance();
6689 if (state_maxsipull->
sinStereo() != 0) {
6690 const auto [covEigenValueSmall, covStereoAngle] = principalComponentAnalysis2x2(covmat);
6691 newerror[0] = std::sqrt(covEigenValueSmall);
6692 newsinstereo = std::sin(covStereoAngle);
6694 newerror[0] = std::sqrt(covmat(0, 0));
6697 const double cosstereo = (newsinstereo == 0) ? 1. : std::sqrt(1 - newsinstereo * newsinstereo);
6699 if (cosstereo != 1.) {
6701 cosstereo * (broadrot->localParameters()[
Trk::locX] - trackpar_maxsipull->parameters()[
Trk::locX]) +
6702 newsinstereo * (broadrot->localParameters()[
Trk::locY] - trackpar_maxsipull->parameters()[
Trk::locY])
6705 newres1 = broadrot->localParameters()[
Trk::locX] - trackpar_maxsipull->parameters()[
Trk::locX];
6708 if (newerror[0] == 0.0) {
6709 ATH_MSG_WARNING(
"Measurement error is zero or negative, treating as outlier");
6712 newpull1 = std::abs(newres1 / newerror[0]);
6716 newerror[1] = std::sqrt(covmat(1, 1));
6717 newres2 = broadrot->localParameters()[
Trk::locY] - trackpar_maxsipull->parameters()[
Trk::locY];
6718 newpull2 = std::abs(newres2 / newerror[1]);
6721 newpull = std::max(newpull1, newpull2);
6727 (newerror[0] > 1.5 * olderror[0] || newerror[1] > 1.5 * std::abs(olderror[1]))
6729 if ((measno_maxsipull < 0) or(measno_maxsipull >= (
int)
res.size())) {
6730 throw std::runtime_error(
6731 "'res' array index out of range in TrkGlobalChi2Fitter/src/GlobalChi2Fitter.cxx:" + std::to_string(__LINE__)
6736 newtrajectory = oldtrajectory;
6738 if (
a.cols() != nfitpars) {
6742 const double oldres1 =
res[measno_maxsipull];
6743 res[measno_maxsipull] = newres1;
6744 err[measno_maxsipull] = newerror[0];
6746 for (
int i = 0; i < nfitpars; i++) {
6747 if (weightderiv(measno_maxsipull, i) == 0) {
6751 b[i] -= weightderiv(measno_maxsipull, i) * (oldres1 / olderror[0] - (newres1 * olderror[0]) / (newerror[0] * newerror[0]));
6753 for (
int j = i; j < nfitpars; j++) {
6757 weightderiv(measno_maxsipull, i) *
6758 weightderiv(measno_maxsipull, j) *
6759 ((olderror[0] * olderror[0]) / (newerror[0] * newerror[0]) - 1)
6763 weightderiv(measno_maxsipull, i) *= olderror[0] / newerror[0];
6767 const double oldres2 =
res[measno_maxsipull + 1];
6768 res[measno_maxsipull + 1] = newres2;
6769 err[measno_maxsipull + 1] = newerror[1];
6771 for (
int i = 0; i < nfitpars; i++) {
6772 if (weightderiv(measno_maxsipull + 1, i) == 0) {
6776 b[i] -= weightderiv(measno_maxsipull + 1, i) * (oldres2 / olderror[1] - (newres2 * olderror[1]) / (newerror[1] * newerror[1]));
6778 for (
int j = i; j < nfitpars; j++) {
6782 weightderiv(measno_maxsipull + 1, i) *
6783 weightderiv(measno_maxsipull + 1, j) *
6784 ((olderror[1] * olderror[1]) / (newerror[1] * newerror[1]) - 1)
6789 weightderiv(measno_maxsipull + 1, i) *= olderror[1] / newerror[1];
6794 "Recovering outlier, hitno=" << hitno_maxsipull <<
" measno=" <<
6795 measno_maxsipull <<
" pull=" << maxsipull <<
" olderror_0=" <<
6796 olderror[0] <<
" newerror_0=" << newerror[0] <<
" olderror_1=" <<
6797 olderror[1] <<
" newerror_1=" << newerror[1]
6806 ((n3sigma < 2 && maxsipull > cut2 && maxsipull < cut) || n3sigma > 1) &&
6816 "Removing outlier, hitno=" << hitno_maxsipull <<
", measno=" <<
6817 measno_maxsipull <<
" pull=" << maxsipull
6824 cleanup_newtrajectory = std::make_unique<GXFTrajectory>(*oldtrajectory);
6825 newtrajectory = cleanup_newtrajectory.get();
6827 if (newa.cols() != nfitpars) {
6833 if ((measno_maxsipull < 0) or(measno_maxsipull >= (
int)
res.size())) {
6834 throw std::runtime_error(
6835 "'res' array index out of range in TrkGlobalChi2Fitter/src/GlobalChi2Fitter.cxx:" + std::to_string(__LINE__)
6839 const double oldres1 =
res[measno_maxsipull];
6840 newres[measno_maxsipull] = 0;
6842 for (
int i = 0; i < nfitpars; i++) {
6843 if (weightderiv(measno_maxsipull, i) == 0) {
6847 newb[i] -= weightderiv(measno_maxsipull, i) * oldres1 / olderror[0];
6849 for (
int j = i; j < nfitpars; j++) {
6853 weightderiv(measno_maxsipull, i) *
6854 weightderiv(measno_maxsipull, j)
6858 newweightderiv(measno_maxsipull, i) = 0;
6862 const double oldres2 =
res[measno_maxsipull + 1];
6863 newres[measno_maxsipull + 1] = 0;
6865 for (
int i = 0; i < nfitpars; i++) {
6866 if (weightderiv(measno_maxsipull + 1, i) == 0) {
6870 newb[i] -= weightderiv(measno_maxsipull + 1, i) * oldres2 / olderror[1];
6872 for (
int j = i; j < nfitpars; j++) {
6873 if (weightderiv(measno_maxsipull + 1, j) == 0) {
6880 weightderiv(measno_maxsipull + 1, i) *
6881 weightderiv(measno_maxsipull + 1, j)
6885 newweightderiv(measno_maxsipull + 1, i) = 0;
6889 newtrajectory->
setOutlier(stateno_maxsipull);
6903 for (
int it = 0; it <
m_maxit; ++it) {
6913 ctx,
cache, *newtrajectory, it, *newap, *newbp, lu_m, doderiv);
6929 const double oldchi2 = oldtrajectory->
chi2() / oldtrajectory->
nDOF();
6930 const double newchi2 = (newtrajectory->
nDOF() > 0) ? newtrajectory->
chi2() / newtrajectory->
nDOF() : 0;
6933 if (newtrajectory->
nDOF() != oldtrajectory->
nDOF() && maxsipull > cut2) {
6934 mindiff = (oldchi2 > .33 *
m_chi2cut || noutl > 0) ? .8 : 1.;
6936 if (noutl == 0 && maxsipull < cut - .5 && oldchi2 < .5 *
m_chi2cut) {
6941 if (newchi2 > oldchi2 || (newchi2 > oldchi2 - mindiff && newchi2 > .33 * oldchi2)) {
6942 ATH_MSG_DEBUG(
"Outlier not confirmed, keeping old trajectory");
6950 (void)cleanup_oldtrajectory.release();
6951 return oldtrajectory;
6953 if (oldtrajectory != newtrajectory) {
6954 cleanup_oldtrajectory = std::move(cleanup_newtrajectory);
6955 oldtrajectory = newtrajectory;
6962 Eigen::LLT < Eigen::MatrixXd >
const lltOfW(
a);
6963 if (lltOfW.info() == Eigen::Success) {
6967 const int ncols =
a.cols();
6968 Amg::MatrixX const weightInvAMG = Amg::MatrixX::Identity(ncols, ncols);
6969 fullcov = lltOfW.solve(weightInvAMG);
6987 oldtrajectory->
nDOF() > 0 &&
6996 (void)cleanup_oldtrajectory.release();
6997 return oldtrajectory;
7009 if (
const auto *pMeas{hit->measurement()};
7015 nrealmeas += hit->numberOfMeasuredParameters();
7019 cache.m_derivmat.setZero();
7025 if (
const auto *pMeas{hit->measurement()};
7031 for (
int i = measindex; i < measindex + hit->numberOfMeasuredParameters(); i++) {
7033 cache.m_derivmat(i, j) = derivs(measindex2, j) * errors[measindex2];
7034 if ((j == 4 && !oldtrajectory.
m_straightline) || j >= nperpars + 2 * nscat) {
7035 cache.m_derivmat(i, j) *= 1000;
7042 measindex += hit->numberOfMeasuredParameters();
7043 }
else if (hit->materialEffects() ==
nullptr) {
7044 measindex2 += hit->numberOfMeasuredParameters();
7050 const EventContext & ctx,
7058 std::unique_ptr<const TrackParameters> per(
nullptr);
7061 std::unique_ptr<const TrackParameters> prevpar(
7066 std::vector<std::pair<const Layer *, const Layer *>> & upstreamlayers = oldtrajectory.
upstreamMaterialLayers();
7069 for (
const auto & [layer1, layer2] : upstreamlayers | std::views::reverse) {
7070 if (prevpar ==
nullptr) {
7075 const Layer *layer = layer1 !=
nullptr ? layer1 : layer2;
7077 const DistanceSolution distsol = layer->surfaceRepresentation().straightLineDistanceEstimate(
7078 prevpar->position(), prevpar->momentum().unit()
7080 const double distance = getDistance(distsol);
7083 if (std::abs(distance) < 0.01) {
7087 if (distsol.
first() * distsol.
second() < 0 && !first) {
7092 if (first && distance > 0) {
7096 std::unique_ptr<const TrackParameters> layerpar(
7100 layer->surfaceRepresentation(),
7108 if (layerpar ==
nullptr) {
7112 if (layer->surfaceRepresentation().bounds().inside(layerpar->localPosition())) {
7116 prevpar = std::move(layerpar);
7129 if (startfactor > 0.5) {
7130 std::unique_ptr<const TrackParameters> updatedpar =
m_matupdator->update(
7134 if (updatedpar !=
nullptr) {
7154 if (endfactor > 0.5) {
7155 std::unique_ptr<const TrackParameters> updatedpar =
m_matupdator->update(
7159 if (updatedpar !=
nullptr) {
7165 if (prevpar !=
nullptr) {
7177 if (per ==
nullptr) {
7178 ATH_MSG_DEBUG(
"Failed to extrapolate to perigee, returning 0");
7197 std::unique_ptr<GXFTrackState>
7199 const EventContext & ctx,
7206 if (per ==
nullptr) {
7210 ATH_MSG_DEBUG(
"Final perigee: " << *per <<
" pos: " << per->position() <<
" pT: " << per->pT());
7216 const std::vector<std::unique_ptr<TrackParameters>> & hc,
7217 std::set<Identifier> & id_set,
7218 std::set<Identifier> & sct_set,
7228 for (
const std::unique_ptr<TrackParameters> & tp : hc) {
7234 if (tp ==
nullptr) {
7243 const TrkDetElementBase * de = tp->associatedSurface().associatedDetectorElement();
7245 if (de ==
nullptr) {
7256 if (id_set.find(
id) != id_set.end()) {
7275 }
else if (
m_DetID->is_sct(
id)) {
7286 }
else if (
m_DetID->is_sct(
id)) {
7295 const Identifier os = e->otherSide()->identify();
7305 if (sct_set.find(os) != sct_set.end()) {
7306 ++rv.m_sct_double_hole;
7335 for (
const std::unique_ptr<GXFTrackState> & s : trajectory.
trackStates()) {
7347 std::vector<std::reference_wrapper<GXFTrackState>> rv;
7355 for (
const std::unique_ptr<GXFTrackState> & s : trajectory.
trackStates()) {
7369 rv.emplace_back(*s);
7376 const TrkDetElementBase * de = s->trackParameters()->associatedSurface().associatedDetectorElement();
7378 if (de !=
nullptr) {
7391 if (s.get() == lastmeas) {
7401 const EventContext & ctx,
7402 const std::vector<std::reference_wrapper<GXFTrackState>> & states
7414 constexpr uint min_meas = 3;
7415 if (std::count_if(states.begin(), states.end(), [](
const GXFTrackState & s){ return s.getStateType(TrackStateOnSurface::Measurement); }) < min_meas) {
7419 bool seen_meas =
false;
7421 std::set<Identifier> id_set;
7422 std::set<Identifier> sct_set;
7428 for (std::size_t i = 0; i < states.size() - 1; i++) {
7451 const double dist = (beg.trackParameters()->
position() - end.trackParameters()->position()).norm();
7453 const bool zStartValid = std::abs(beg.trackParameters()->position().z())<10000.;
7455 ATH_MSG_DEBUG(
"Pathological track parameter well outside of detector");
7456 ATH_MSG_DEBUG(
"Propagator might have issue with this, skipping");
7457 ATH_MSG_VERBOSE(
"dumping track parameters " << *(beg.trackParameters()));
7466 if (seen_meas && dist >= 2.5 && zStartValid) {
7473 std::optional<std::vector<std::unique_ptr<TrackParameters>>> & hc = beg.getHoles();
7474 std::vector<std::unique_ptr<TrackParameters>> states;
7482 if (hc.has_value()) {
7483 states = std::move(*hc);
7521 std::vector<std::unique_ptr<Trk::TrackParameters>>
const bl =
m_extrapolator->extrapolateBlindly(
7542 const EventContext & ctx,
7548 auto trajectory = std::make_unique<Trk::TrackStates>();
7558 if (perigee_ts ==
nullptr) {
7564 trajectory->reserve(tmptrajectory.
trackStates().size());
7570 hit->resetTrackCovariance();
7576 hit->materialEffects()))
7582 auto trackState = hit->trackStateOnSurface();
7583 hit->resetTrackCovariance();
7584 trajectory->emplace_back(trackState.release());
7587 auto qual = std::make_unique<FitQuality>(tmptrajectory.
chi2(), tmptrajectory.
nDOF());
7607 std::unique_ptr<Track> rv = std::make_unique<Track>(info, std::move(trajectory), std::move(qual));
7617 std::unique_ptr<TrackSummary>
ts = std::make_unique<TrackSummary>();
7626 std::optional<TrackHoleCount> hole_count;
7633 std::vector<std::reference_wrapper<GXFTrackState>>
const states =
holeSearchStates(tmptrajectory);
7651 if (hole_count.has_value()) {
7664 rv->setTrackSummary(std::move(
ts));
7673 const EventContext & ctx,
7682 std::vector<std::unique_ptr<TrackParameters>> rv =
m_extrapolator->extrapolateStepwise(
7696 &rv.front()->associatedSurface() == &src.associatedSurface() ||
7697 trackParametersClose(*rv.front(), src, 0.001) ||
7701 rv.front().reset(
nullptr);
7712 &rv.back()->associatedSurface() == &src.associatedSurface() ||
7713 trackParametersClose(*rv.back(), src, 0.001) ||
7717 rv.back().reset(
nullptr);
7724 const EventContext & ctx,
7732 std::unique_ptr<const TrackParameters> rv;
7733 std::optional<TransportJacobian> jac{};
7744 if (rv !=
nullptr && calcderiv) {
7749 std::optional<std::vector<std::unique_ptr<TrackParameters>>> extrapolation;
7758 std::move(extrapolation)
7763 const EventContext & ctx,
7774 ctx, prev,
ts, propdir, bf, calcderiv, holesearch
7777 if (rv.m_parameters ==
nullptr) {
7778 propdir = invertPropdir(propdir);
7781 ctx, prev,
ts, propdir, bf, calcderiv, holesearch
7789 const EventContext& ctx,
7796 std::vector<std::unique_ptr<GXFTrackState>> & states = trajectory.
trackStates();
7799 std::unique_ptr<const TrackParameters> tmptrackpar;
7801 for (
int hitno = nstatesupstream - 1; hitno >= 0; hitno--) {
7802 const Surface &surf1 = states[hitno]->associatedSurface();
7809 const double distance = getDistance(distsol);
7831 (rv.m_parameters !=
nullptr) &&
7832 (prevtrackpar->
position() - rv.m_parameters->position()).mag() > 5 * mm
7838 if (rv.m_parameters ==
nullptr) {
7839 ATH_MSG_DEBUG(
"propagation failed, prev par: " << *prevtrackpar <<
7840 " pos: " << prevtrackpar->
position() <<
" destination surface: " << surf1);
7844 states[hitno]->setTrackParameters(std::move(rv.m_parameters));
7845 const TrackParameters *currenttrackpar = states[hitno]->trackParameters();
7846 const Surface &surf = states[hitno]->associatedSurface();
7848 if (rv.m_jacobian != std::nullopt) {
7850 states[hitno]->materialEffects() !=
nullptr &&
7851 states[hitno]->materialEffects()->deltaE() != 0 &&
7852 states[hitno]->materialEffects()->sigmaDeltaE() <= 0 &&
7855 const double p = 1. / std::abs(currenttrackpar->parameters()[
Trk::qOverP]);
7856 const double de = std::abs(states[hitno]->materialEffects()->deltaE());
7857 const double mass = trajectory.
mass();
7858 const double newp = std::sqrt(p * p + 2 * de * std::sqrt(mass * mass + p * p) + de * de);
7859 (*rv.m_jacobian) (4, 4) = ((p + p * de / std::sqrt(p * p + mass * mass)) / newp) * p * p / (newp * newp);
7862 states[hitno]->setJacobian(*rv.m_jacobian);
7863 }
else if (calcderiv) {
7870 if (meff !=
nullptr && hitno != 0) {
7872 surf, *meff, *states[hitno]->trackParameters(), trajectory.
mass(), -1
7875 if (std::holds_alternative<FitterStatusCode>(
r)) {
7876 return std::get<FitterStatusCode>(
r);
7879 tmptrackpar = std::move(std::get<std::unique_ptr<const TrackParameters>>(
r));
7880 prevtrackpar = tmptrackpar.get();
7882 prevtrackpar = currenttrackpar;
7888 for (
int hitno = nstatesupstream; hitno < (int) states.size(); hitno++) {
7889 const Surface &surf = states[hitno]->associatedSurface();
7893 const double distance = getDistance(distsol);
7910 (rv.m_parameters !=
nullptr) &&
7912 (prevtrackpar->
position() - rv.m_parameters->position()).mag() > 5 * mm
7917 if (rv.m_parameters ==
nullptr) {
7918 ATH_MSG_DEBUG(
"propagation failed, prev par: " << *prevtrackpar <<
7919 " pos: " << prevtrackpar->
7920 position() <<
" destination surface: " << surf);
7924 if (rv.m_jacobian != std::nullopt) {
7926 states[hitno]->materialEffects() !=
nullptr &&
7927 states[hitno]->materialEffects()->deltaE() != 0 &&
7928 states[hitno]->materialEffects()->sigmaDeltaE() <= 0 &&
7931 const double p = 1 / std::abs(rv.m_parameters->parameters()[
Trk::qOverP]);
7932 const double de = std::abs(states[hitno]->materialEffects()->deltaE());
7933 const double mass = trajectory.
mass();
7934 double newp = p * p - 2 * de * std::sqrt(mass * mass + p * p) + de * de;
7937 newp = std::sqrt(newp);
7940 (*rv.m_jacobian) (4, 4) = ((p - p * de / std::sqrt(p * p + mass * mass)) / newp) * p * p / (newp * newp);
7943 states[hitno]->setJacobian(*rv.m_jacobian);
7944 }
else if (calcderiv) {
7951 if (meff !=
nullptr) {
7953 surf, *meff, *rv.m_parameters, trajectory.
mass(), +1
7956 if (std::holds_alternative<FitterStatusCode>(
r)) {
7957 return std::get<FitterStatusCode>(
r);
7960 rv.m_parameters = std::move(std::get<std::unique_ptr<const TrackParameters>>(
r));
7963 states[hitno]->setTrackParameters(std::move(rv.m_parameters));
7964 prevtrackpar = states[hitno]->trackParameters();
7977 const AmgVector(5) & old = param.parameters();
7982 if (!correctAngles(newphi, newtheta)) {
7983 ATH_MSG_DEBUG(
"Angles out of range, phi: " << newphi <<
" theta: " << newtheta);
7987 double newqoverp = 0;
7993 const double oldp = std::abs(1 / old[
Trk::qOverP]);
7994 const double newp2 = oldp * oldp -
sign * 2 * std::abs(meff.
deltaE()) * std::sqrt(mass * mass + oldp * oldp) + meff.
deltaE() * meff.
deltaE();
8001 newqoverp = std::copysign(1 / std::sqrt(newp2), old[
Trk::qOverP]);
8008 old[0], old[1], newphi, newtheta, newqoverp, std::nullopt
8020 using Matrix55 = Eigen::Matrix<double, 5, 5>;
8022 Matrix55 initialjac;
8023 initialjac.setZero();
8024 initialjac(4, 4) = 1;
8026 Matrix55 jacvertex(initialjac);
8028 std::vector<Matrix55, Eigen::aligned_allocator<Matrix55>> jacscat(trajectory.
numberOfScatterers(), initialjac);
8029 std::vector<Matrix55, Eigen::aligned_allocator<Matrix55>> jacbrem(trajectory.
numberOfBrems(), initialjac);
8031 std::vector<std::unique_ptr<GXFTrackState>> & states = trajectory.
trackStates();
8040 for (
const bool forward : {
false,
true}) {
8042 hit_begin = nstatesupstream;
8043 hit_end = (int) states.size();
8044 scatno = nscatupstream;
8045 bremno = nbremupstream;
8047 hit_begin = nstatesupstream - 1;
8054 int hitno = hit_begin;
8055 forward ? (hitno < hit_end) : (hitno >= hit_end);
8056 hitno += (forward ? 1 : -1)
8059 state = states[hitno].get();
8063 if (fillderivmat && state->
derivatives().cols() != nfitpars) {
8071 const int jmaxbrem = 4;
8073 if (hitno == (forward ? hit_end - 1 : 0)) {
8074 if (!fillderivmat) {
8082 Eigen::Matrix<double, 5, 5> & jac = state->
jacobian();
8084 if (hitno == nstatesupstream + (forward ? 0 : -1)) {
8085 jacvertex.block<4, 5>(0, 0) = jac.block<4, 5>(0, 0);
8086 jacvertex(4, 4) = jac(4, 4);
8096 jcnt = jmax - jmin + 1;
8098 lp_bgn = forward ? nscatupstream : nscatupstream - 1;
8101 for (
int i = lp_bgn; forward ? (i < lp_end) : (i > lp_end); i += (forward ? 1 : -1)) {
8103 i == scatno + (forward ? -1 : 1) &&
8104 prevstate !=
nullptr &&
8108 jacscat[i].block(0, jmin, 4, jcnt) = jac.block(0, jmin, 4, jcnt);
8109 jacscat[i](4, 4) = jac(4, 4);
8111 calculateJac(jac, jacscat[i], jmin, jmax);
8115 Eigen::MatrixXd & derivmat = state->
derivatives();
8116 const int scatterPos = nperpars + 2 * i;
8118 derivmat.block<4, 2>(0, scatterPos) = (forward ? 1 : -1) * jacscat[i].block<4, 2>(0, 2);
8124 jcnt = jmax - jmin + 1;
8126 lp_bgn = forward ? nbremupstream : nbremupstream - 1;
8129 for (
int i = lp_bgn; forward ? (i < lp_end) : (i > lp_end); i += (forward ? 1 : -1)) {
8131 i == bremno + (forward ? -1 : 1) &&
8136 jacbrem[i].block(0, jmin, 4, jcnt) = jac.block(0, jmin, 4, jcnt);
8137 jacbrem[i](4, 4) = jac(4, 4);
8139 calculateJac(jac, jacbrem[i], jmin, jmax);
8143 Eigen::MatrixXd & derivmat = state->
derivatives();
8144 const int scatterPos = nperpars + 2 * nscats + i;
8146 derivmat.block<5, 1>(0, scatterPos) = (forward ? .001 : -.001) * jacbrem[i].block<5, 1>(0, 4);
8150 calculateJac(jac, jacvertex, 0, 4);
8154 Eigen::MatrixXd & derivmat = state->
derivatives();
8155 derivmat.block(0, 0, 4, nperpars) = jacvertex.block(0, 0, 4, nperpars);
8157 if (nperpars == 5) {
8158 derivmat.col(4).segment(0, 4) *= .001;
8159 derivmat(4, 4) = .001 * jacvertex(4, 4);
8165 (!trajectory.
prefit() || states[hitno]->materialEffects()->deltaE() == 0)
8167 scatno += (forward ? 1 : -1);
8171 states[hitno]->materialEffects() &&
8172 states[hitno]->materialEffects()->sigmaDeltaE() > 0
8174 bremno += (forward ? 1 : -1);
8177 prevstate = states[hitno].get();
8186 bool onlylocal)
const {
8192 std::vector<std::unique_ptr<GXFTrackState>> & states = trajectory.
trackStates();
8194 std::vector < int >
indices(states.size());
8196 int i = nstatesupstream;
8197 for (
int j = 0; j < (int) states.size(); j++) {
8198 if (j < nstatesupstream) {
8205 for (
int stateno = 0; stateno < (int) states.size(); stateno++) {
8206 if (stateno == 0 || stateno == nstatesupstream) {
8207 prevstate =
nullptr;
8210 std::unique_ptr<GXFTrackState> & state = states[
index];
8211 if (state->materialEffects() !=
nullptr) {
8212 prevstate = state.get();
8216 if (!state->hasTrackCovariance()) {
8217 state->zeroTrackCovariance();
8219 AmgMatrix(5, 5) & trackerrmat = state->trackCovariance();
8221 if ((prevstate !=
nullptr) &&
8225 Eigen::Matrix<double, 5, 5> & jac = state->jacobian();
8226 const AmgMatrix(5, 5)& prevcov = states[
indices[stateno - 1]]->trackCovariance();
8228 trackerrmat = jac * prevcov * jac.transpose();
8232 trackerrmat = derivatives * fullcovmat * derivatives.transpose();
8242 bool errorok =
true;
8243 for (
int i = 0; i < 5; i++) {
8246 && trackerrmat(i, i) > meascov(j, j)) {
8248 const double scale = std::sqrt(meascov(j, j) / trackerrmat(i, i));
8249 trackerrmat(i, i) = meascov(j, j);
8250 for (
int k = 0; k < 5; k++) {
8252 trackerrmat(k, i) *= scale;
8260 for (
int i = 0; i < 5; i++) {
8264 for (
int j = 0; j < 5; j++) {
8272 trackerrmat(4, 4) = 1e-20;
8276 state->trackParameters();
8278 std::optional<
AmgMatrix(5, 5)> trkerrmat;
8280 if (state->hasTrackCovariance()) {
8281 trkerrmat = (state->trackCovariance());
8283 trkerrmat = std::nullopt;
8286 const AmgVector(5) & tpars = tmptrackpar->parameters();
8287 std::unique_ptr<const TrackParameters> trackpar(
8293 std::move(trkerrmat))
8295 state->setTrackParameters(std::move(trackpar));
8298 if (errorok && trajectory.
nDOF() > 0) {
8299 fitQual =
m_updator->fullStateFitQuality(
8300 *state->trackParameters(),
8308 state->setFitQuality(fitQual);
8310 prevstate = state.get();
8314 std::optional<TransportJacobian>
8316 const EventContext& ctx,
8329 std::optional<TransportJacobian> jac = std::make_optional<TransportJacobian>(J);
8332 0.01, 0.01, 0.00001, 0.00001, 0.000000001
8343 for (
int i = 0; i < 5; i++) {
8346 if (thisdiscsurf && i == 1) {
8352 if (i == 0 && thiscylsurf) {
8353 vecminuseps[i] = -std::remainder(-vecminuseps[i], 2 *
M_PI * previousSurface.
bounds().
r());
8354 }
else if (i == 1 && thisdiscsurf) {
8355 vecpluseps[i] = -std::remainder(-vecpluseps[i], 2 *
M_PI);
8360 std::unique_ptr<const TrackParameters> parpluseps(
8370 const std::unique_ptr<const TrackParameters> parminuseps(
8381 std::unique_ptr<const TrackParameters> newparpluseps(
8392 std::unique_ptr<const TrackParameters> newparminuseps(
8407 if (newparpluseps ==
nullptr) {
8419 if (newparminuseps ==
nullptr) {
8431 if ((newparpluseps ==
nullptr) || (newparminuseps ==
nullptr)) {
8435 for (
int j = 0; j < 5; j++) {
8439 if (j == 0 && cylsurf) {
8441 }
else if (j == 1 && discsurf) {
8445 (*jac) (j, i) =
diff / (2 * eps[i]);
8455 auto nmeas1 = pDataVector->size();
8456 const auto *pLastValue = (*pDataVector)[nmeas1 - 1];
8464 testrot =
static_cast<const RIO_OnTrack *
>(pLastValue);
8466 if (lastMeasIsCompetingRIO){
8468 testrot = &testcrot->rioOnTrack(0);
8472 if (testrot ==
nullptr) {
8473 const auto *pPenultimate = (*pDataVector)[nmeas1 - 2];
8476 if(penultimateIsRIO){
8477 testrot =
static_cast<const RIO_OnTrack *
>(pPenultimate);
8479 if (penultimateIsCompetingRIO){
8481 testrot = &testcrot->rioOnTrack(0);
8488 (testrot !=
nullptr) &&
8505 if (cond_obj ==
nullptr) {
8514 std::stringstream
msg;
8516 throw std::runtime_error(
msg.str());
8520 if (
cache.m_caloEntrance ==
nullptr) {
8523 if (geometry !=
nullptr) {
8524 cache.m_caloEntrance = geometry->trackingVolume(
"InDet::Containers::InnerDetector");
8532 if (
cache.m_caloEntrance ==
nullptr) {
8537 return cache.m_caloEntrance !=
nullptr;
8541 if (
cache.m_msEntrance ==
nullptr) {
8544 if (geometry !=
nullptr) {
8545 cache.m_msEntrance = geometry->trackingVolume(
"MuonSpectrometerEntrance");
8553 if (
cache.m_msEntrance ==
nullptr) {
8558 return cache.m_msEntrance !=
nullptr;
Scalar perp() const
perp method - perpendicular length
Scalar deltaPhi(const MatrixBase< Derived > &vec) const
Scalar phi() const
phi method
Scalar theta() const
theta method
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
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
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
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.