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 = std::move(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;
825 const std::unique_ptr<const TrackParameters> origlastidpar =
unique_clone(lastidpar);
827 double newqoverpid = 0;
829 if (tp_closestmuon and !firstismuon) {
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 = std::move(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];
5156 auto *pMattEff = state->materialEffects();
5167 doNewPseudoMeasurements &&
5169 !state->associatedSurface().isFree() &&
5170 !state->isRecalibrated()
5173 covMatrix(0, 0) = 100;
5175 std::unique_ptr<const PseudoMeasurementOnTrack> newpseudo = std::make_unique<const PseudoMeasurementOnTrack>(
5177 std::move(covMatrix),
5181 state->setMeasurement(std::move(newpseudo));
5182 measbase = state->measurement();
5189 double *errors = state->measurementErrors();
5191 for (
int i = 0; i < 5; i++) {
5211 res[measno] = residuals[i];
5217 res[measno] = -std::remainder(-
res[measno], 2 *
M_PI);
5226 double *errors = state->measurementErrors();
5227 for (
int i = 0; i < 5; i++) {
5228 if (errors[i] > 0) {
5229 error[measno] = errors[i];
5240 ((trajectory.
prefit() == 0) || pMattEff->deltaE() == 0)
5242 const double deltaPhi = pMattEff->deltaPhi();
5243 const double measDeltaPhi = pMattEff->measuredDeltaPhi();
5244 const double sigma2deltaPhi = std::pow(pMattEff->sigmaDeltaPhi(), 2);
5245 const double deltaTheta = pMattEff->deltaTheta();
5246 const double sigma2deltaTheta = std::pow(pMattEff->sigmaDeltaTheta(), 2);
5248 if (trajectory.
prefit() != 1) {
5249 b[nperpars + 2 * scatno] -= (
deltaPhi - measDeltaPhi) / sigma2deltaPhi;
5250 b[nperpars + 2 * scatno + 1] -= deltaTheta / sigma2deltaTheta;
5252 b[nperpars + scatno] -= deltaTheta / sigma2deltaTheta;
5257 deltaTheta * deltaTheta / sigma2deltaTheta
5266 if ((pMattEff !=
nullptr) && pMattEff->sigmaDeltaE() > 0) {
5267 double averagenergyloss = std::abs(pMattEff->deltaE());
5268 const double qoverpbrem = limitInversePValue(1000 * states[hitno]->trackParameters()->parameters()[
Trk::qOverP]);
5269 const double qoverp = limitInversePValue(qoverpbrem - pMattEff->delta_p());
5270 const double pbrem = 1. / std::abs(qoverpbrem);
5271 const double p = 1. / std::abs(qoverp);
5272 const double mass = .001 * trajectory.
mass();
5273 const double energy = std::sqrt(p * p + mass * mass);
5274 const double bremEnergy = std::sqrt(pbrem * pbrem + mass * mass);
5276 const double resMaterial = .001 * averagenergyloss - energy + bremEnergy;
5277 res[nmeas - nbrem + bremno] = resMaterial;
5279 const double sigde = pMattEff->sigmaDeltaE();
5280 const double sigdepos = pMattEff->sigmaDeltaEPos();
5281 const double sigdeneg = pMattEff->sigmaDeltaENeg();
5283 double errorMaterial = .001 * pMattEff->sigmaDeltaE();
5284 error[nmeas - nbrem + bremno] = errorMaterial;
5295 if (pMattEff->isKink()) {
5296 maxbrempull = -999999999;
5297 state_maxbrempull =
nullptr;
5301 cache.m_asymeloss &&
5303 trajectory.
prefit() == 0 &&
5305 sigde != sigdepos &&
5308 const double elosspull = resMaterial / errorMaterial;
5310 if (trajectory.
mass() > 100) {
5315 if (std::abs(elosspull) > 1) {
5316 if (elosspull < -1) {
5317 pMattEff->setSigmaDeltaE(sigdepos);
5319 pMattEff->setSigmaDeltaE(sigdeneg);
5322 errorMaterial = .001 * pMattEff->sigmaDeltaE();
5323 error[nmeas - nbrem + bremno] = errorMaterial;
5335 !pMattEff->isKink() && (
5336 (
m_fixbrem == -1 && elosspull < maxbrempull) ||
5340 bremno_maxbrempull = bremno;
5341 state_maxbrempull = state.get();
5342 maxbrempull = elosspull;
5351 trajectory.
prefit() == 0 &&
5352 pMattEff->sigmaDeltaPhi() == 0 &&
5353 pMattEff->isMeasuredEloss() &&
5354 resMaterial / (.001 * pMattEff->sigmaDeltaEAve()) > 2.5
5356 const TrackParameters* parforcalo = states[hitno - 2]->trackParameters();
5359 std::vector<MaterialEffectsOnTrack> calomeots =
5372 if (calomeots.size() == 3) {
5373 averagenergyloss = std::abs(calomeots[1].energyLoss()->deltaE());
5374 const double newres = .001 * averagenergyloss - energy + bremEnergy;
5375 const double newerr = .001 * calomeots[1].energyLoss()->sigmaDeltaE();
5377 const double oldPull = resMaterial / errorMaterial;
5378 const double newPull = newres / newerr;
5380 if (std::abs(newPull) < std::abs(oldPull)) {
5381 ATH_MSG_DEBUG(
"Changing from measured to parametrized energy loss");
5383 pMattEff->setEloss(std::unique_ptr<EnergyLoss>(calomeots[1].energyLoss()->clone()));
5384 pMattEff->setSigmaDeltaE(calomeots[1].energyLoss()->sigmaDeltaE());
5385 res[nmeas - nbrem + bremno] = newres;
5386 error[nmeas - nbrem + bremno] = newerr;
5390 pMattEff->setMeasuredEloss(
false);
5400 for (
int imeas = 0; imeas < nmeas; imeas++) {
5401 if (
error[imeas] == 0) {
5422 const double oldChi2 = trajectory.
prevchi2();
5423 const double newChi2 = trajectory.
chi2();
5428 const double nDOF = trajectory.
nDOF();
5429 const double oldRedChi2 = (nDOF > 0) ? oldChi2 / nDOF : 0;
5430 const double newRedChi2 = (nDOF > 0) ? newChi2 / nDOF : 0;
5433 trajectory.
prefit() > 0 && (
5434 (newRedChi2 < 2 && it != 0) ||
5435 (newRedChi2 < oldRedChi2 + .1 && std::abs(newRedChi2 - oldRedChi2) < 1 && it != 1)
5448 int miniter = (nsihits != 0 && nsihits + ntrthits == nhits) ? 1 : 2;
5449 miniter = std::max(miniter,
cache.m_miniter);
5451 if (it >= miniter && std::abs(oldChi2 - newChi2) < 1) {
5458 const int bremno_maxbrempull,
5464 if (state_maxbrempull ==
nullptr) {
5476 const int nmeas = (int)
res.size();
5479 const double oldError =
error[nmeas - nbrem + bremno_maxbrempull];
5481 error[nmeas - nbrem + bremno_maxbrempull] = newError;
5484 if (
a.cols() != nFitPars) {
5488 const double errorRatio = oldError / newError;
5489 const double errorReductionRatio = 1 - std::pow(errorRatio, 2);
5492 for (
int i = 0; i < nFitPars; i++) {
5493 if (weightderiv(nmeas - nbrem + bremno_maxbrempull, i) == 0) {
5497 for (
int j = i; j < nFitPars; j++) {
5498 const double newaij =
a(i, j) - errorReductionRatio *
5499 weightderiv(nmeas - nbrem + bremno_maxbrempull, i) *
5500 weightderiv(nmeas - nbrem + bremno_maxbrempull, j);
5502 a.fillSymmetric(i, j, newaij);
5504 weightderiv(nmeas - nbrem + bremno_maxbrempull, i) *= errorRatio;
5513 std::vector<std::unique_ptr<GXFTrackState>> & states = trajectory.
trackStates();
5526 const int nmeas = (int) weightderiv.rows();
5528 for (std::unique_ptr<GXFTrackState> & state : states) {
5532 const auto [scatmin, scatmax] = std::minmax(scatno, nscatupstream);
5533 const auto [bremmin, bremmax] = std::minmax(bremno, nbremupstream);
5540 const double sinStereo =
5542 state->sinStereo() :
5544 const double cosStereo =
5546 std::sqrt(1 - std::pow(sinStereo, 2)) :
5554 auto getThisDeriv = [sinStereo, cosStereo, &derivatives](
int i,
int j) ->
double {
5555 if (i == 0 && sinStereo != 0) {
5556 return derivatives(0, j) * cosStereo + sinStereo * derivatives(1, j);
5558 return derivatives(i, j);
5562 for (
int i = 0; i < 5; i++) {
5577 if (i == 0 && sinStereo != 0) {
5578 weightderiv.row(measno).head(cols) =
5579 (derivatives.row(0).
head(cols) * cosStereo +
5580 sinStereo * derivatives.row(1).head(cols)) /
5583 weightderiv.row(measno).head(cols) = derivatives.row(i).head(cols) /
error[measno];
5587 for (
int j = scatmin; j < scatmax; j++) {
5588 if (trajectory.
prefit() == 1) {
5589 const int index = nperparams + j;
5592 const int index = nperparams + 2 * j;
5594 weightderiv(measno,
index + 1) = getThisDeriv(i,
index + 1) /
error[measno];
5598 for (
int j = bremmin; j < bremmax; j++) {
5599 const int index = j + nperparams + 2 * nscat;
5606 double *errors = state->measurementErrors();
5607 for (
int i = 0; i < 5; i++) {
5608 if (errors[i] > 0) {
5614 ((trajectory.
prefit() == 0) || state->materialEffects()->deltaE() == 0)
5619 if ((state->materialEffects() !=
nullptr) && state->materialEffects()->sigmaDeltaE() > 0) {
5621 const double qoverpbrem = limitInversePValue(1000 * state->trackParameters()->parameters()[
Trk::qOverP]);
5622 const double qoverp = limitInversePValue(qoverpbrem - state->materialEffects()->delta_p());
5624 const double mass = .001 * trajectory.
mass();
5626 const auto thisMeasurementIdx{nmeas - nbrem + bremno};
5630 auto multiplier = [] (
double mass,
double qOverP){
5633 const auto qoverpTerm {multiplier(mass, qoverp) /
error[thisMeasurementIdx]};
5634 const auto qoverpBremTerm {multiplier(mass, qoverpbrem) /
error[thisMeasurementIdx]};
5637 weightderiv(thisMeasurementIdx, 4) = qoverpBremTerm - qoverpTerm;
5640 const auto bremNoBase = nperparams + 2 * nscat;
5641 if (bremno < nbremupstream) {
5642 weightderiv(thisMeasurementIdx, bremNoBase + bremno) = qoverpTerm;
5643 for (
int bremno2 = bremno + 1; bremno2 < nbremupstream; bremno2++) {
5644 weightderiv(thisMeasurementIdx, bremNoBase + bremno2) = qoverpTerm - qoverpBremTerm;
5647 weightderiv(thisMeasurementIdx, bremNoBase + bremno) = qoverpBremTerm;
5648 for (
int bremno2 = nbremupstream; bremno2 < bremno; bremno2++) {
5649 weightderiv(thisMeasurementIdx, bremNoBase + bremno2) = qoverpBremTerm - qoverpTerm;
5668 const int nMeas = (int)
res.size();
5670 cache.m_firstmeasurement.resize(nFitPars);
5671 cache.m_lastmeasurement.resize(nFitPars);
5673 for (
int i = 0; i < nPerPars; i++) {
5674 cache.m_firstmeasurement[i] = 0;
5675 cache.m_lastmeasurement[i] = nMeas - nBrem;
5682 const std::unique_ptr<GXFTrackState> & state = trajectory.
trackStates()[i];
5685 if (meff ==
nullptr) {
5686 measno += state->numberOfMeasuredParameters();
5690 const int firstMeasurement = i < nUpstreamStates ? 0 : measno;
5691 const int lastMeasurement = i < nUpstreamStates ? measno : nMeas - nBrem;
5694 && (trajectory.
prefit() == 0 || meff->
deltaE() == 0)) {
5695 const int scatterPos = nPerPars + 2 * scatno;
5697 cache.m_firstmeasurement[scatterPos] = firstMeasurement;
5698 cache.m_lastmeasurement[scatterPos] = lastMeasurement;
5700 cache.m_firstmeasurement[scatterPos + 1] = firstMeasurement;
5701 cache.m_lastmeasurement[scatterPos + 1] = lastMeasurement;
5707 const int bremPos = nPerPars + nScatPars + bremno;
5709 cache.m_firstmeasurement[bremPos] = firstMeasurement;
5710 cache.m_lastmeasurement[bremPos] = lastMeasurement;
5731 const int nMeas = (int)
res.size();
5733 for (
int k = 0; k < nFitPars; k++) {
5734 const int minMeasK =
cache.m_firstmeasurement[k];
5735 const int maxMeasK =
cache.m_lastmeasurement[k];
5742 for (
int measno = minMeasK; measno < maxMeasK; measno++) {
5743 b[k] +=
res[measno] * (1. /
error[measno]) * weightDeriv(measno, k);
5753 if (k == 4 || k >= nPerPars + nScatPars) {
5754 for (
int measno = nMeas - nBrem; measno < nMeas; measno++) {
5755 b[k] +=
res[measno] * (1. /
error[measno]) * weightDeriv(measno, k);
5769 for (
int k = 0; k < nFitPars; k++) {
5770 for (
int l = k; l < nFitPars; l++) {
5771 const int minMeas = std::max(
cache.m_firstmeasurement[k],
cache.m_firstmeasurement[l]);
5772 const int maxMeas = std::min(
cache.m_lastmeasurement[k],
cache.m_lastmeasurement[l]);
5775 for (
int measno = minMeas; measno < maxMeas; measno++) {
5776 a_kl += weightDeriv(measno, k) * weightDeriv(measno, l);
5779 a.fillSymmetric(l, k, a_kl);
5797 const int nMeas = (int)
res.size();
5804 for (
int k = nPerPars; k < nPerPars + nScatPars; k += 2) {
5805 a(k, k) += 1. / std::pow(scatSigmas[scatno].first, 2);
5806 a(k + 1, k + 1) += 1. / std::pow(scatSigmas[scatno].second, 2);
5814 for (
int measno = nMeas - nBrem; measno < nMeas; measno++) {
5815 for (
int k = 4; k < nFitPars; k++) {
5817 k = nPerPars + nScatPars;
5820 for (
int l = k; l < nFitPars; l++) {
5822 l = nPerPars + nScatPars;
5825 const double a_kl =
a(l, k) + weightDeriv(measno, k) * weightDeriv(measno, l);
5826 a.fillSymmetric(l, k, a_kl);
5838 const double oldRedChi2,
5839 const double newRedChi2
5847 bool weightChanged =
false;
5854 double newPhiWeight = 1.1;
5855 double newThetaWeight = 1.001;
5856 if (trajectory.
prefit() == 0) {
5862 newPhiWeight = 1.00000001;
5863 }
else if (it == 1) {
5864 newPhiWeight = 1.0000001;
5865 }
else if (it <= 3) {
5866 newPhiWeight = 1.0001;
5867 }
else if (it <= 6) {
5868 newPhiWeight = 1.01;
5871 if (newRedChi2 > oldRedChi2 - 1 && newRedChi2 < oldRedChi2) {
5872 newPhiWeight = 1.0001;
5873 newThetaWeight = 1.0001;
5874 }
else if (newRedChi2 > oldRedChi2 - 25 && newRedChi2 < oldRedChi2) {
5875 newPhiWeight = 1.001;
5876 newThetaWeight = 1.0001;
5883 std::size_t scatno = 0;
5888 for (
const auto & state : trajectory.
trackStates()) {
5891 if (meff ==
nullptr) {
5895 const bool isValidPlaneSurface =
5897 static_cast<const PlaneSurface *
>(&state->associatedSurface()) !=
nullptr;
5902 if (meff->
deltaE() == 0 || (trajectory.
prefit() == 0 && isValidPlaneSurface)) {
5903 weightChanged =
true;
5905 const int scatNoIndex = 2 * scatno + nPerPars;
5908 if (scatno >=
cache.m_phiweight.size()) {
5909 std::stringstream message;
5910 message <<
"scatno is out of range " << scatno <<
" !< " <<
cache.m_phiweight.size();
5911 throw std::range_error(message.str());
5919 a(scatNoIndex, scatNoIndex) /=
cache.m_phiweight[scatno];
5922 cache.m_phiweight[scatno] = newPhiWeight;
5923 a(scatNoIndex, scatNoIndex) *= newPhiWeight;
5924 }
else if (trajectory.
prefit() >= 2) {
5925 a(scatNoIndex, scatNoIndex) *= newPhiWeight;
5926 a(scatNoIndex + 1, scatNoIndex + 1) *= newThetaWeight;
5953 trajectory.
prefit() == 2 &&
5956 (newRedChi2 < oldRedChi2 - 25 || newRedChi2 > oldRedChi2)
5961 return weightChanged;
5970 std::size_t scatno = 0;
5979 if (scatno >=
cache.m_phiweight.size()) {
5980 std::stringstream message;
5981 message <<
"scatno is out of range " << scatno <<
" !< " <<
cache.m_phiweight.size();
5982 throw std::range_error(message.str());
5985 const bool isValidPlaneSurface =
5987 static_cast<const PlaneSurface *
>(&state->associatedSurface()) !=
nullptr;
5989 if (meff->
deltaE() == 0 || isValidPlaneSurface) {
5990 const int scatNoIndex = 2 * scatno + nPerPars;
5991 a(scatNoIndex, scatNoIndex) /=
cache.m_phiweight[scatno];
5992 cache.m_phiweight[scatno] = 1;
6006 const EventContext& ctx,
6015 const int nDOFold = trajectory.
nDOF();
6016 const double oldChi2 = trajectory.
chi2();
6017 const double oldRedChi2 = nDOFold > 0 ? oldChi2 / nDOFold : 0;
6019 if (
cache.m_phiweight.empty()) {
6039 int bremno_maxbrempull = 0;
6055 if ((state_maxbrempull !=
nullptr) && trajectory.
converged()) {
6064 const int nDOFnew = trajectory.
nDOF();
6065 const double newChi2 = trajectory.
chi2();
6066 const double newRedChi2 = nDOFnew > 0 ? newChi2 / nDOFnew : 0;
6068 ATH_MSG_DEBUG(
"old chi2: " << oldChi2 <<
"/" << nDOFold <<
"=" << oldRedChi2 <<
6069 ", new chi2: " << newChi2 <<
"/" << nDOFnew <<
"=" << newRedChi2);
6080 if (
cache.m_firstmeasurement.empty()) {
6104 if (doDeriv || weightChanged) {
6115 if (trajectory.
prefit() == 0) {
6121 if (nSiHits + nTrtHits !=
nHits) {
6128 (newRedChi2 < 2 || (newRedChi2 < oldRedChi2 && newRedChi2 > oldRedChi2 - .5))
6150 Eigen::LLT<Eigen::MatrixXd>
const llt(lu_m);
6152 if (llt.info() != Eigen::Success) {
6176 double d0 = refpar->parameters()[
Trk::d0];
6177 double z0 = refpar->parameters()[
Trk::z0];
6180 double qoverp = refpar->parameters()[
Trk::qOverP];
6182 if (nperparams > 0) {
6183 d0 += deltaParameters[0];
6184 z0 += deltaParameters[1];
6185 phi += deltaParameters[2];
6186 theta += deltaParameters[3];
6187 qoverp = (trajectory.
m_straightline) ? 0 : .001 * deltaParameters[4] + qoverp;
6199 std::vector < std::pair < double, double >>&scatangles = trajectory.
scatteringAngles();
6200 for (
int i = 0; i < nscat; i++) {
6201 scatangles[i].first += deltaParameters[2 * i + nperparams];
6202 scatangles[i].second += deltaParameters[2 * i + nperparams + 1];
6208 std::vector < double >&delta_ps = trajectory.
brems();
6209 for (
int i = 0; i < nbrem; i++) {
6210 delta_ps[i] += deltaParameters[nperparams + 2 * nscat + i];
6216 std::unique_ptr<const TrackParameters> newper(
6236 const EventContext& evtctx
6247 if (!splitProbContainer.
isValid()) {
6251 std::vector<std::unique_ptr<GXFTrackState>> & states = trajectory.
trackStates();
6258 for (
size_t stateno = 0; stateno < states.size(); stateno++) {
6263 measno += states[stateno-1]->numberOfMeasuredParameters();
6266 std::unique_ptr<GXFTrackState> & state = states[stateno];
6278 const auto *
const rot =
static_cast<const RIO_OnTrack *
>(pMeas);
6279 prd = rot->prepRawData();
6289 const auto &splitProb = splitProbContainer->splitProbability(
pixelCluster);
6290 if (!splitProb.isSplit()) {
6291 ATH_MSG_DEBUG(
"Pixel cluster is not split so no need to update" );
6295 std::unique_ptr < const RIO_OnTrack > newrot;
6296 double *olderror = state->measurementErrors();
6299 double newerror[5] = {-1,-1,-1,-1,-1};
6300 double newres[2] = {-1,-1};
6302 newrot.reset(
m_ROTcreator->correct(*prd, *trackpars, evtctx));
6307 const Amg::MatrixX & covmat = newrot->localCovariance();
6309 newerror[0] = std::sqrt(covmat(0, 0));
6310 newres[0] = newrot->localParameters()[
Trk::locX] - trackpars->parameters()[
Trk::locX];
6311 newerror[1] = std::sqrt(covmat(1, 1));
6312 newres[1] = newrot->localParameters()[
Trk::locY] - trackpars->parameters()[
Trk::locY];
6314 if (
a.cols() != nfitpars) {
6319 for(
int k =0; k<2; k++ ){
6320 const double oldres =
res[measno+k];
6321 res[measno+k] = newres[k];
6322 err[measno+k] = newerror[k];
6324 for (
int i = 0; i < nfitpars; i++) {
6325 if (weightderiv(measno+k, i) == 0) {
6329 b[i] -= weightderiv(measno+k, i) * (oldres / olderror[k] - (newres[k] * olderror[k]) / (newerror[k] * newerror[k]));
6331 for (
int j = i; j < nfitpars; j++) {
6335 weightderiv(measno+k, i) *
6336 weightderiv(measno+k, j) *
6337 ((olderror[k] * olderror[k]) / (newerror[k] * newerror[k]) - 1)
6341 weightderiv(measno+k, i) *= olderror[k] / newerror[k];
6345 state->setMeasurement(std::move(newrot));
6346 state->setMeasurementErrors(newerror);
6361 const EventContext& ctx
6369 std::vector<std::unique_ptr<GXFTrackState>> & states = trajectory.
trackStates();
6375 if (
a.cols() != nfitpars) {
6383 bool outlierremoved =
false;
6384 bool hitrecalibrated =
false;
6386 for (
int stateno = 0; stateno < (int) states.size(); stateno++) {
6387 std::unique_ptr<GXFTrackState> & state = states[stateno];
6395 std::abs(state->trackParameters()->parameters()[
Trk::driftRadius]) > 1.05 * state->associatedSurface().bounds().r()
6400 outlierremoved =
true;
6402 double *errors = state->measurementErrors();
6403 const double olderror = errors[0];
6407 for (
int i = 0; i < nfitpars; i++) {
6408 if (weightderiv(measno, i) == 0) {
6412 b[i] -=
res[measno] * weightderiv(measno, i) / olderror;
6414 for (
int j = i; j < nfitpars; j++) {
6417 a(i, j) - weightderiv(measno, i) * weightderiv(measno, j)
6420 weightderiv(measno, i) = 0;
6424 }
else if (trtrecal) {
6425 double *errors = state->measurementErrors();
6426 const double olderror = errors[0];
6428 const auto *
const thisMeasurement{state->measurement()};
6437 const double trackradius = state->trackParameters()->parameters()[
Trk::driftRadius];
6439 std::unique_ptr<const Trk::RIO_OnTrack> newrot =
nullptr;
6440 const double distance = std::abs(std::abs(trackradius) - dcradius);
6442 if (distance < scalefactor * dcerror && (olderror > 1. || trackradius * oldradius < 0)) {
6444 }
else if (distance > scalefactor * dcerror && olderror < 1.) {
6448 if (newrot !=
nullptr) {
6450 hitrecalibrated =
true;
6454 if ((measno < 0) or (measno >= (
int)
res.size())) {
6455 throw std::runtime_error(
6456 "'res' array index out of range in TrkGlobalChi2Fitter/src/GlobalChi2Fitter.cxx:" + std::to_string(__LINE__)
6460 const double oldres =
res[measno];
6461 const double newres = newradius - state->trackParameters()->parameters()[
Trk::driftRadius];
6462 errors[0] = newerror;
6463 state->setMeasurement(std::move(newrot));
6467 for (
int i = 0; i < nfitpars; i++) {
6468 if (weightderiv(measno, i) == 0) {
6472 b[i] -= weightderiv(measno, i) * (oldres / olderror - (newres * olderror) / (newerror * newerror));
6474 for (
int j = i; j < nfitpars; j++) {
6478 !
cache.m_phiweight.empty() &&
6481 i < nperpars + 2 * nscats &&
6482 (i - nperpars) % 2 == 0
6484 weight =
cache.m_phiweight[(i - nperpars) / 2];
6489 a(i, j) + weightderiv(measno, i) * weightderiv(measno, j) * ((olderror * olderror) / (newerror * newerror) - 1) * weight
6492 weightderiv(measno, i) *= olderror / newerror;
6495 res[measno] = newres;
6496 err[measno] = newerror;
6505 measno += state->numberOfMeasuredParameters();
6509 if (trajectory.
nDOF() < 0) {
6514 if (outlierremoved || hitrecalibrated) {
6518 cache.m_miniter = it + 2;
6523 const EventContext& ctx,
6531 bool trackok =
false;
6533 std::unique_ptr < GXFTrajectory > cleanup_oldtrajectory;
6535 std::unique_ptr < GXFTrajectory > cleanup_newtrajectory;
6541 while (!trackok && oldtrajectory->
nDOF() > 0) {
6543 std::vector<std::unique_ptr<GXFTrackState>> & states = oldtrajectory->
trackStates();
6551 if (nhits != nsihits) {
6555 double maxsipull = -1;
6557 int hitno_maxsipull = -1;
6558 int measno_maxsipull = -1;
6559 int stateno_maxsipull = 0;
6571 for (
int stateno = 0; stateno < (int) states.size(); stateno++) {
6572 std::unique_ptr<GXFTrackState> & state = states[stateno];
6578 double *errors = state->measurementErrors();
6580 const Amg::MatrixX & hitcov = state->measurement()->localCovariance();
6581 const double sinstereo = state->sinStereo();
6582 const double cosstereo = (sinstereo == 0) ? 1 : std::sqrt(1 - sinstereo * sinstereo);
6583 double weight1 = -1;
6585 if (hitcov(0, 0) > trackcov(0, 0)) {
6586 if (sinstereo == 0) {
6587 weight1 = errors[0] * errors[0] - trackcov(0, 0);
6589 weight1 = errors[0] * errors[0] - (
6590 trackcov(0, 0) * cosstereo * cosstereo + 2 *
6591 trackcov(1, 0) * cosstereo * sinstereo + trackcov(1, 1) * sinstereo * sinstereo
6596 const double weight2 = (
6598 errors[1] * errors[1] - trackcov(1, 1) :
6602 double sipull1 = weight1 > 0 ? std::abs(
res[measno] / std::sqrt(weight1)) : -1;
6603 const double sipull2 = (
6605 std::abs(
res[measno + 1] / std::sqrt(weight2)) :
6608 sipull1 = std::max(sipull1, sipull2);
6610 if (sipull1 > maxsipull) {
6611 maxsipull = sipull1;
6612 measno_maxsipull = measno;
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;
6956 a = std::move(newa);
6957 b = std::move(newb);
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 &&
7002 (void)cleanup_oldtrajectory.release();
7003 return oldtrajectory;
7015 if (
const auto *pMeas{
hit->measurement()};
7021 nrealmeas +=
hit->numberOfMeasuredParameters();
7025 cache.m_derivmat.setZero();
7031 if (
const auto *pMeas{
hit->measurement()};
7037 for (
int i = measindex; i < measindex +
hit->numberOfMeasuredParameters(); i++) {
7039 cache.m_derivmat(i, j) = derivs(measindex2, j) * errors[measindex2];
7040 if ((j == 4 && !oldtrajectory.
m_straightline) || j >= nperpars + 2 * nscat) {
7041 cache.m_derivmat(i, j) *= 1000;
7048 measindex +=
hit->numberOfMeasuredParameters();
7049 }
else if (
hit->materialEffects() ==
nullptr) {
7050 measindex2 +=
hit->numberOfMeasuredParameters();
7056 const EventContext & ctx,
7064 std::unique_ptr<const TrackParameters> per(
nullptr);
7067 std::unique_ptr<const TrackParameters> prevpar(
7072 std::vector<std::pair<const Layer *, const Layer *>> & upstreamlayers = oldtrajectory.
upstreamMaterialLayers();
7075 for (
const auto & [layer1, layer2] : upstreamlayers | std::views::reverse) {
7076 if (prevpar ==
nullptr) {
7081 const Layer *layer = layer1 !=
nullptr ? layer1 : layer2;
7083 const DistanceSolution distsol = layer->surfaceRepresentation().straightLineDistanceEstimate(
7084 prevpar->position(), prevpar->momentum().unit()
7086 const double distance = getDistance(distsol);
7089 if (std::abs(distance) < 0.01) {
7093 if (distsol.
first() * distsol.
second() < 0 && !first) {
7098 if (first && distance > 0) {
7102 std::unique_ptr<const TrackParameters> layerpar(
7106 layer->surfaceRepresentation(),
7114 if (layerpar ==
nullptr) {
7118 if (layer->surfaceRepresentation().bounds().inside(layerpar->localPosition())) {
7122 prevpar = std::move(layerpar);
7135 if (startfactor > 0.5) {
7136 std::unique_ptr<const TrackParameters> updatedpar =
m_matupdator->update(
7140 if (updatedpar !=
nullptr) {
7160 if (endfactor > 0.5) {
7161 std::unique_ptr<const TrackParameters> updatedpar =
m_matupdator->update(
7165 if (updatedpar !=
nullptr) {
7171 if (prevpar !=
nullptr) {
7183 if (per ==
nullptr) {
7184 ATH_MSG_DEBUG(
"Failed to extrapolate to perigee, returning 0");
7203 std::unique_ptr<GXFTrackState>
7205 const EventContext & ctx,
7212 if (per ==
nullptr) {
7216 ATH_MSG_DEBUG(
"Final perigee: " << *per <<
" pos: " << per->position() <<
" pT: " << per->pT());
7222 const std::vector<std::unique_ptr<TrackParameters>> & hc,
7223 std::set<Identifier> & id_set,
7224 std::set<Identifier> & sct_set,
7234 for (
const std::unique_ptr<TrackParameters> & tp : hc) {
7240 if (tp ==
nullptr) {
7249 const TrkDetElementBase * de = tp->associatedSurface().associatedDetectorElement();
7251 if (de ==
nullptr) {
7262 if (id_set.find(
id) != id_set.end()) {
7281 }
else if (
m_DetID->is_sct(
id)) {
7292 }
else if (
m_DetID->is_sct(
id)) {
7301 const Identifier os = e->otherSide()->identify();
7311 if (sct_set.find(os) != sct_set.end()) {
7312 ++rv.m_sct_double_hole;
7341 for (
const std::unique_ptr<GXFTrackState> & s : trajectory.
trackStates()) {
7353 std::vector<std::reference_wrapper<GXFTrackState>> rv;
7361 for (
const std::unique_ptr<GXFTrackState> & s : trajectory.
trackStates()) {
7375 rv.emplace_back(*s);
7382 const TrkDetElementBase * de = s->trackParameters()->associatedSurface().associatedDetectorElement();
7384 if (de !=
nullptr) {
7397 if (s.get() == lastmeas) {
7407 const EventContext & ctx,
7408 const std::vector<std::reference_wrapper<GXFTrackState>> & states
7420 constexpr uint min_meas = 3;
7421 if (std::count_if(states.begin(), states.end(), [](
const GXFTrackState & s){ return s.getStateType(TrackStateOnSurface::Measurement); }) < min_meas) {
7425 bool seen_meas =
false;
7427 std::set<Identifier> id_set;
7428 std::set<Identifier> sct_set;
7434 for (std::size_t i = 0; i < states.size() - 1; i++) {
7457 const double dist = (beg.trackParameters()->
position() - end.trackParameters()->position()).norm();
7459 const bool zStartValid = std::abs(beg.trackParameters()->position().z())<10000.;
7461 ATH_MSG_DEBUG(
"Pathological track parameter well outside of detector");
7462 ATH_MSG_DEBUG(
"Propagator might have issue with this, skipping");
7463 ATH_MSG_VERBOSE(
"dumping track parameters " << *(beg.trackParameters()));
7472 if (seen_meas && dist >= 2.5 && zStartValid) {
7479 std::optional<std::vector<std::unique_ptr<TrackParameters>>> & hc = beg.getHoles();
7480 std::vector<std::unique_ptr<TrackParameters>> theseStates;
7488 if (hc.has_value()) {
7489 theseStates = std::move(*hc);
7527 std::vector<std::unique_ptr<Trk::TrackParameters>>
const bl =
m_extrapolator->extrapolateBlindly(
7548 const EventContext & ctx,
7554 auto trajectory = std::make_unique<Trk::TrackStates>();
7564 if (perigee_ts ==
nullptr) {
7570 trajectory->reserve(tmptrajectory.
trackStates().size());
7576 hit->resetTrackCovariance();
7582 hit->materialEffects()))
7588 auto trackState =
hit->trackStateOnSurface();
7589 hit->resetTrackCovariance();
7590 trajectory->emplace_back(trackState.release());
7593 auto qual = std::make_unique<FitQuality>(tmptrajectory.
chi2(), tmptrajectory.
nDOF());
7613 std::unique_ptr<Track> rv = std::make_unique<Track>(info, std::move(trajectory), std::move(qual));
7623 std::unique_ptr<TrackSummary>
ts = std::make_unique<TrackSummary>();
7632 std::optional<TrackHoleCount> hole_count;
7639 std::vector<std::reference_wrapper<GXFTrackState>>
const states =
holeSearchStates(tmptrajectory);
7657 if (hole_count.has_value()) {
7670 rv->setTrackSummary(std::move(
ts));
7679 const EventContext & ctx,
7688 std::vector<std::unique_ptr<TrackParameters>> rv =
m_extrapolator->extrapolateStepwise(
7702 &rv.front()->associatedSurface() == &src.associatedSurface() ||
7703 trackParametersClose(*rv.front(), src, 0.001) ||
7707 rv.front().reset(
nullptr);
7718 &rv.back()->associatedSurface() == &src.associatedSurface() ||
7719 trackParametersClose(*rv.back(), src, 0.001) ||
7723 rv.back().reset(
nullptr);
7730 const EventContext & ctx,
7738 std::unique_ptr<const TrackParameters> rv;
7739 std::optional<TransportJacobian> jac{};
7750 if (rv !=
nullptr && calcderiv) {
7755 std::optional<std::vector<std::unique_ptr<TrackParameters>>> extrapolation;
7764 std::move(extrapolation)
7769 const EventContext & ctx,
7780 ctx, prev,
ts, propdir, bf, calcderiv, holesearch
7783 if (rv.m_parameters ==
nullptr) {
7784 propdir = invertPropdir(propdir);
7787 ctx, prev,
ts, propdir, bf, calcderiv, holesearch
7795 const EventContext& ctx,
7802 std::vector<std::unique_ptr<GXFTrackState>> & states = trajectory.
trackStates();
7805 std::unique_ptr<const TrackParameters> tmptrackpar;
7807 for (
int hitno = nstatesupstream - 1; hitno >= 0; hitno--) {
7808 const Surface &surf1 = states[hitno]->associatedSurface();
7815 const double distance = getDistance(distsol);
7837 (rv.m_parameters !=
nullptr) &&
7838 (prevtrackpar->
position() - rv.m_parameters->position()).mag() > 5 * mm
7844 if (rv.m_parameters ==
nullptr) {
7845 ATH_MSG_DEBUG(
"propagation failed, prev par: " << *prevtrackpar <<
7846 " pos: " << prevtrackpar->
position() <<
" destination surface: " << surf1);
7850 states[hitno]->setTrackParameters(std::move(rv.m_parameters));
7851 const TrackParameters *currenttrackpar = states[hitno]->trackParameters();
7852 const Surface &surf = states[hitno]->associatedSurface();
7854 if (rv.m_jacobian != std::nullopt) {
7856 states[hitno]->materialEffects() !=
nullptr &&
7857 states[hitno]->materialEffects()->deltaE() != 0 &&
7858 states[hitno]->materialEffects()->sigmaDeltaE() <= 0 &&
7861 const double p = 1. / std::abs(currenttrackpar->parameters()[
Trk::qOverP]);
7862 const double de = std::abs(states[hitno]->materialEffects()->deltaE());
7863 const double mass = trajectory.
mass();
7864 const double newp = std::sqrt(p * p + 2 * de * std::sqrt(mass * mass + p * p) + de * de);
7865 (*rv.m_jacobian) (4, 4) = ((p + p * de / std::sqrt(p * p + mass * mass)) / newp) * p * p / (newp * newp);
7868 states[hitno]->setJacobian(*rv.m_jacobian);
7869 }
else if (calcderiv) {
7876 if (meff !=
nullptr && hitno != 0) {
7878 surf, *meff, *states[hitno]->trackParameters(), trajectory.
mass(), -1
7881 if (std::holds_alternative<FitterStatusCode>(
r)) {
7882 return std::get<FitterStatusCode>(
r);
7885 tmptrackpar = std::move(std::get<std::unique_ptr<const TrackParameters>>(
r));
7886 prevtrackpar = tmptrackpar.get();
7888 prevtrackpar = currenttrackpar;
7894 for (
int hitno = nstatesupstream; hitno < (int) states.size(); hitno++) {
7895 const Surface &surf = states[hitno]->associatedSurface();
7899 const double distance = getDistance(distsol);
7916 (rv.m_parameters !=
nullptr) &&
7918 (prevtrackpar->
position() - rv.m_parameters->position()).mag() > 5 * mm
7923 if (rv.m_parameters ==
nullptr) {
7924 ATH_MSG_DEBUG(
"propagation failed, prev par: " << *prevtrackpar <<
7925 " pos: " << prevtrackpar->
7926 position() <<
" destination surface: " << surf);
7930 if (rv.m_jacobian != std::nullopt) {
7932 states[hitno]->materialEffects() !=
nullptr &&
7933 states[hitno]->materialEffects()->deltaE() != 0 &&
7934 states[hitno]->materialEffects()->sigmaDeltaE() <= 0 &&
7937 const double p = 1 / std::abs(rv.m_parameters->parameters()[
Trk::qOverP]);
7938 const double de = std::abs(states[hitno]->materialEffects()->deltaE());
7939 const double mass = trajectory.
mass();
7940 double newp = p * p - 2 * de * std::sqrt(mass * mass + p * p) + de * de;
7943 newp = std::sqrt(newp);
7946 (*rv.m_jacobian) (4, 4) = ((p - p * de / std::sqrt(p * p + mass * mass)) / newp) * p * p / (newp * newp);
7949 states[hitno]->setJacobian(*rv.m_jacobian);
7950 }
else if (calcderiv) {
7957 if (meff !=
nullptr) {
7959 surf, *meff, *rv.m_parameters, trajectory.
mass(), +1
7962 if (std::holds_alternative<FitterStatusCode>(
r)) {
7963 return std::get<FitterStatusCode>(
r);
7966 rv.m_parameters = std::move(std::get<std::unique_ptr<const TrackParameters>>(
r));
7969 states[hitno]->setTrackParameters(std::move(rv.m_parameters));
7970 prevtrackpar = states[hitno]->trackParameters();
7983 const AmgVector(5) & old = param.parameters();
7988 if (!correctAngles(newphi, newtheta)) {
7989 ATH_MSG_DEBUG(
"Angles out of range, phi: " << newphi <<
" theta: " << newtheta);
7993 double newqoverp = 0;
7999 const double oldp = std::abs(1 / old[
Trk::qOverP]);
8000 const double newp2 = oldp * oldp -
sign * 2 * std::abs(meff.
deltaE()) * std::sqrt(mass * mass + oldp * oldp) + meff.
deltaE() * meff.
deltaE();
8007 newqoverp = std::copysign(1 / std::sqrt(newp2), old[
Trk::qOverP]);
8014 old[0], old[1], newphi, newtheta, newqoverp, std::nullopt
8026 using Matrix55 = Eigen::Matrix<double, 5, 5>;
8028 Matrix55 initialjac;
8029 initialjac.setZero();
8030 initialjac(4, 4) = 1;
8032 Matrix55 jacvertex(initialjac);
8034 std::vector<Matrix55, Eigen::aligned_allocator<Matrix55>> jacscat(trajectory.
numberOfScatterers(), initialjac);
8035 std::vector<Matrix55, Eigen::aligned_allocator<Matrix55>> jacbrem(trajectory.
numberOfBrems(), initialjac);
8037 std::vector<std::unique_ptr<GXFTrackState>> & states = trajectory.
trackStates();
8046 for (
const bool forward : {
false,
true}) {
8048 hit_begin = nstatesupstream;
8049 hit_end = (int) states.size();
8050 scatno = nscatupstream;
8051 bremno = nbremupstream;
8053 hit_begin = nstatesupstream - 1;
8060 int hitno = hit_begin;
8061 forward ? (hitno < hit_end) : (hitno >= hit_end);
8062 hitno += (forward ? 1 : -1)
8065 state = states[hitno].get();
8069 if (fillderivmat && state->
derivatives().cols() != nfitpars) {
8077 const int jmaxbrem = 4;
8079 if (hitno == (forward ? hit_end - 1 : 0)) {
8080 if (!fillderivmat) {
8088 Eigen::Matrix<double, 5, 5> & jac = state->
jacobian();
8090 if (hitno == nstatesupstream + (forward ? 0 : -1)) {
8091 jacvertex.block<4, 5>(0, 0) = jac.block<4, 5>(0, 0);
8092 jacvertex(4, 4) = jac(4, 4);
8102 jcnt = jmax - jmin + 1;
8104 lp_bgn = forward ? nscatupstream : nscatupstream - 1;
8107 for (
int i = lp_bgn; forward ? (i < lp_end) : (i > lp_end); i += (forward ? 1 : -1)) {
8109 i == scatno + (forward ? -1 : 1) &&
8110 prevstate !=
nullptr &&
8114 jacscat[i].block(0, jmin, 4, jcnt) = jac.block(0, jmin, 4, jcnt);
8115 jacscat[i](4, 4) = jac(4, 4);
8117 calculateJac(jac, jacscat[i], jmin, jmax);
8121 Eigen::MatrixXd & derivmat = state->
derivatives();
8122 const int scatterPos = nperpars + 2 * i;
8124 derivmat.block<4, 2>(0, scatterPos) = (forward ? 1 : -1) * jacscat[i].block<4, 2>(0, 2);
8130 jcnt = jmax - jmin + 1;
8132 lp_bgn = forward ? nbremupstream : nbremupstream - 1;
8135 for (
int i = lp_bgn; forward ? (i < lp_end) : (i > lp_end); i += (forward ? 1 : -1)) {
8137 i == bremno + (forward ? -1 : 1) &&
8142 jacbrem[i].block(0, jmin, 4, jcnt) = jac.block(0, jmin, 4, jcnt);
8143 jacbrem[i](4, 4) = jac(4, 4);
8145 calculateJac(jac, jacbrem[i], jmin, jmax);
8149 Eigen::MatrixXd & derivmat = state->
derivatives();
8150 const int scatterPos = nperpars + 2 * nscats + i;
8152 derivmat.block<5, 1>(0, scatterPos) = (forward ? .001 : -.001) * jacbrem[i].block<5, 1>(0, 4);
8156 calculateJac(jac, jacvertex, 0, 4);
8160 Eigen::MatrixXd & derivmat = state->
derivatives();
8161 derivmat.block(0, 0, 4, nperpars) = jacvertex.block(0, 0, 4, nperpars);
8163 if (nperpars == 5) {
8164 derivmat.col(4).segment(0, 4) *= .001;
8165 derivmat(4, 4) = .001 * jacvertex(4, 4);
8171 (!trajectory.
prefit() || states[hitno]->materialEffects()->deltaE() == 0)
8173 scatno += (forward ? 1 : -1);
8177 states[hitno]->materialEffects() &&
8178 states[hitno]->materialEffects()->sigmaDeltaE() > 0
8180 bremno += (forward ? 1 : -1);
8183 prevstate = states[hitno].get();
8192 bool onlylocal)
const {
8198 std::vector<std::unique_ptr<GXFTrackState>> & states = trajectory.
trackStates();
8200 std::vector < int >
indices(states.size());
8202 int i = nstatesupstream;
8203 for (
int j = 0; j < (int) states.size(); j++) {
8204 if (j < nstatesupstream) {
8211 for (
int stateno = 0; stateno < (int) states.size(); stateno++) {
8212 if (stateno == 0 || stateno == nstatesupstream) {
8213 prevstate =
nullptr;
8216 std::unique_ptr<GXFTrackState> & state = states[
index];
8217 if (state->materialEffects() !=
nullptr) {
8218 prevstate = state.get();
8222 if (!state->hasTrackCovariance()) {
8223 state->zeroTrackCovariance();
8225 AmgMatrix(5, 5) & trackerrmat = state->trackCovariance();
8227 if ((prevstate !=
nullptr) &&
8231 Eigen::Matrix<double, 5, 5> & jac = state->jacobian();
8232 const AmgMatrix(5, 5)& prevcov = states[
indices[stateno - 1]]->trackCovariance();
8234 trackerrmat = jac * prevcov * jac.transpose();
8238 trackerrmat = derivatives * fullcovmat * derivatives.transpose();
8248 bool errorok =
true;
8249 for (
int i = 0; i < 5; i++) {
8252 && trackerrmat(i, i) > meascov(j, j)) {
8254 const double scale = std::sqrt(meascov(j, j) / trackerrmat(i, i));
8255 trackerrmat(i, i) = meascov(j, j);
8256 for (
int k = 0; k < 5; k++) {
8258 trackerrmat(k, i) *= scale;
8266 for (
int i = 0; i < 5; i++) {
8270 for (
int j = 0; j < 5; j++) {
8278 trackerrmat(4, 4) = 1e-20;
8282 state->trackParameters();
8284 std::optional<
AmgMatrix(5, 5)> trkerrmat;
8286 if (state->hasTrackCovariance()) {
8287 trkerrmat = (state->trackCovariance());
8289 trkerrmat = std::nullopt;
8292 const AmgVector(5) & tpars = tmptrackpar->parameters();
8293 std::unique_ptr<const TrackParameters> trackpar(
8299 std::move(trkerrmat))
8301 state->setTrackParameters(std::move(trackpar));
8304 if (errorok && trajectory.
nDOF() > 0) {
8305 fitQual =
m_updator->fullStateFitQuality(
8306 *state->trackParameters(),
8314 state->setFitQuality(fitQual);
8316 prevstate = state.get();
8320 std::optional<TransportJacobian>
8322 const EventContext& ctx,
8335 std::optional<TransportJacobian> jac = std::make_optional<TransportJacobian>(J);
8338 0.01, 0.01, 0.00001, 0.00001, 0.000000001
8349 for (
int i = 0; i < 5; i++) {
8352 if (thisdiscsurf && i == 1) {
8358 if (i == 0 && thiscylsurf) {
8359 vecminuseps[i] = -std::remainder(-vecminuseps[i], 2 *
M_PI * previousSurface.
bounds().
r());
8360 }
else if (i == 1 && thisdiscsurf) {
8361 vecpluseps[i] = -std::remainder(-vecpluseps[i], 2 *
M_PI);
8366 std::unique_ptr<const TrackParameters> parpluseps(
8376 const std::unique_ptr<const TrackParameters> parminuseps(
8387 std::unique_ptr<const TrackParameters> newparpluseps(
8398 std::unique_ptr<const TrackParameters> newparminuseps(
8413 if (newparpluseps ==
nullptr) {
8425 if (newparminuseps ==
nullptr) {
8437 if ((newparpluseps ==
nullptr) || (newparminuseps ==
nullptr)) {
8441 for (
int j = 0; j < 5; j++) {
8445 if (j == 0 && cylsurf) {
8447 }
else if (j == 1 && discsurf) {
8451 (*jac) (j, i) =
diff / (2 * eps[i]);
8461 auto nmeas1 = pDataVector->size();
8462 const auto *pLastValue = (*pDataVector)[nmeas1 - 1];
8470 testrot =
static_cast<const RIO_OnTrack *
>(pLastValue);
8472 if (lastMeasIsCompetingRIO){
8474 testrot = &testcrot->rioOnTrack(0);
8478 if (testrot ==
nullptr) {
8479 const auto *pPenultimate = (*pDataVector)[nmeas1 - 2];
8482 if(penultimateIsRIO){
8483 testrot =
static_cast<const RIO_OnTrack *
>(pPenultimate);
8485 if (penultimateIsCompetingRIO){
8487 testrot = &testcrot->rioOnTrack(0);
8494 (testrot !=
nullptr) &&
8511 if (cond_obj ==
nullptr) {
8520 std::stringstream
msg;
8522 throw std::runtime_error(
msg.str());
8526 if (
cache.m_caloEntrance ==
nullptr) {
8529 if (geometry !=
nullptr) {
8530 cache.m_caloEntrance = geometry->trackingVolume(
"InDet::Containers::InnerDetector");
8538 if (
cache.m_caloEntrance ==
nullptr) {
8543 return cache.m_caloEntrance !=
nullptr;
8547 if (
cache.m_msEntrance ==
nullptr) {
8550 if (geometry !=
nullptr) {
8551 cache.m_msEntrance = geometry->trackingVolume(
"MuonSpectrometerEntrance");
8559 if (
cache.m_msEntrance ==
nullptr) {
8564 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_DEBUG(x,...)
#define ATH_MSG_ERROR(x,...)
#define ATH_MSG_WARNING(x,...)
#define ATH_MSG_VERBOSE(x,...)
#define ATH_MSG_INFO(x,...)
#define ATH_MSG_FATAL(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.
double pT() const
Access method for transverse momentum.
virtual const Surface & associatedSurface() const =0
Access to the Surface associated to the Parameters.
Class describing the Line to which the Perigee refers to.
virtual Surface::ChargedTrackParametersUniquePtr createUniqueTrackParameters(double l1, double l2, double phi, double theta, double qop, const std::optional< AmgSymMatrix(5)> &cov=std::nullopt) const override final
Use the Surface as a ParametersBase constructor, from local parameters - charged.
virtual DistanceSolution straightLineDistanceEstimate(const Amg::Vector3D &pos, const Amg::Vector3D &dir) const override final
fast straight line distance evaluation to Surface
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, const 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.