57 #include "CLHEP/Units/SystemOfUnits.h"
69 #include <Eigen/Dense>
70 #include <Eigen/StdVector>
81 return distsol.
first();
84 std::abs(distsol.
first()) < std::abs(distsol.
second()) ?
95 limitInversePValue(
double qOverP){
96 const double magnitude = std::abs(
qOverP);
98 constexpr
double maxP{100.*10e6*
MeV};
99 constexpr
double minP{1.e-3*
MeV};
100 constexpr
double lo {1./maxP};
101 constexpr
double hi {1./minP};
102 double limited = std::clamp(magnitude, lo, hi);
103 return std::copysign(limited,
qOverP);
107 std::pair<const Trk::TrackParameters *, const Trk::TrackParameters *> getFirstLastIdPar(
const Trk::Track &
track) {
113 while ((firstidpar ==
nullptr) && parit !=
track.trackParameters()->end()) {
115 ((**parit).covariance() !=
nullptr) &&
124 parit =
track.trackParameters()->end();
128 ((**parit).covariance() !=
nullptr) &&
133 }
while ((lastidpar ==
nullptr) && parit !=
track.trackParameters()->begin());
135 return std::make_pair(firstidpar, lastidpar);
152 std::abs(
a.parameters()[0] -
b.parameters()[0]) <
e &&
153 std::abs(
a.parameters()[1] -
b.parameters()[1]) <
e &&
154 std::abs(
a.parameters()[2] -
b.parameters()[2]) <
e
160 coercePhiCoordinateRange(
double &
phi){
161 const auto absPhi{std::abs(
phi)};
162 constexpr
double twoPi{2.*
M_PI};
163 if (std::abs(
phi - twoPi) < absPhi) {
166 if (std::abs(
phi + twoPi) < absPhi) {
178 calculateJac(Eigen::Matrix<double, 5, 5> &jac,
179 Eigen::Matrix<double, 5, 5> &
out,
180 int jmin,
int jmax) {
184 out.block(0, 0, 4, jmin).setZero();
188 out.block(0, jmax + 1, 4, 5 - (jmax + 1)).setZero();
191 out(4, 4) = jac(4, 4);
198 const std::string &
t,
199 const std::string &
n,
227 ATH_MSG_ERROR(
"Hole search requested but no boundary check tool provided.");
228 return StatusCode::FAILURE;
253 ATH_MSG_WARNING(
"FillDerivativeMatrix option selected, switching off acceleration!");
277 ATH_MSG_ERROR(
"Hole search requested but track summaries are disabled.");
278 return StatusCode::FAILURE;
283 return StatusCode::SUCCESS;
289 if (m_fit_status[
S_FITS] > 0) {
292 <<
" track fits failed because of a matrix inversion failure");
294 <<
" tracks were rejected by the outlier logic");
296 <<
" track fits failed because of a propagation failure");
298 <<
" track fits failed because of an invalid angle (theta/phi)");
300 <<
" track fits failed because the fit did not converge");
302 <<
" tracks did not pass the chi^2 cut");
304 <<
" tracks were killed by the energy loss update");
307 return StatusCode::SUCCESS;
312 std::unique_ptr<Track>
314 const EventContext& ctx,
320 ATH_MSG_DEBUG(
"--> entering GlobalChi2Fitter::fit(Track,Track,)");
335 const Track *indettrack = firstismuon ? &intrk2 : &intrk1;
336 const Track *muontrack = firstismuon ? &intrk1 : &intrk2;
338 bool measphi =
false;
345 rot = &crot->rioOnTrack(0);
355 double dotprod2 = measdir.dot(
358 if (std::abs(dotprod1) < .5 && std::abs(dotprod2) < .5) {
366 auto [firstidpar, lastidpar] = getFirstLastIdPar(*indettrack);
368 if ((firstidpar ==
nullptr) || (lastidpar ==
nullptr)) {
372 std::unique_ptr<const TrackParameters> parforcalo =
unique_clone(firstismuon ? firstidpar : lastidpar);
375 const AmgVector(5) & newpars = parforcalo->parameters();
377 parforcalo=parforcalo->associatedSurface().createUniqueTrackParameters(
378 newpars[0], newpars[1], newpars[2], newpars[3], 1 / 5000., std::nullopt
382 std::vector < MaterialEffectsOnTrack > calomeots;
385 calomeots =
m_calotool->extrapolationSurfacesAndEffects(
389 parforcalo->associatedSurface(),
402 if (calomeots.empty()) {
407 std::unique_ptr<Track>
track;
417 double qoverpid = measperid !=
nullptr ? measperid->parameters()[
Trk::qOverP] : 0;
418 double qoverpmuon = measpermuon !=
nullptr ? measpermuon->parameters()[
Trk::qOverP] : 0;
420 const AmgSymMatrix(5) * errmatid = measperid !=
nullptr ? measperid->covariance() :
nullptr;
421 const AmgSymMatrix(5) * errmatmuon = measpermuon !=
nullptr ? measpermuon->covariance() :
nullptr;
425 (errmatid !=
nullptr) &&
426 (errmatmuon !=
nullptr) &&
431 double piderror = std::sqrt((*errmatid) (4, 4)) / (qoverpid * qoverpid);
432 double pmuonerror = std::sqrt((*errmatmuon) (4, 4)) / (qoverpmuon * qoverpmuon);
433 double energyerror = std::sqrt(
434 calomeots[1].energyLoss()->sigmaDeltaE() *
435 calomeots[1].energyLoss()->sigmaDeltaE() + piderror * piderror +
436 pmuonerror * pmuonerror
440 (std::abs(calomeots[1].energyLoss()->deltaE()) -
441 std::abs(1 / qoverpid) + std::abs(1 / qoverpmuon)
444 ATH_MSG_DEBUG(
"Changing from measured to parametrized energy loss");
449 parforcalo->associatedSurface(),
454 if (calomeots.empty()) {
462 bool firstfitwasattempted =
false;
484 qoverpid * qoverpmuon > 0
490 firstfitwasattempted =
true;
495 (
track ==
nullptr) &&
496 !firstfitwasattempted &&
503 trajectory = trajectory2;
507 bool pseudoupdated =
false;
509 if (
track !=
nullptr) {
510 for (std::unique_ptr<GXFTrackState> & pseudostate : trajectory.
trackStates()) {
511 if (pseudostate ==
nullptr) {
522 if ((pseudostate ==
nullptr) || pseudostate->fitQuality().chiSquared() < 10) {
527 std::unique_ptr<const TrackParameters> updpar(
m_updator->removeFromState(
529 pseudostate->measurement()->localParameters(),
530 pseudostate->measurement()->localCovariance()
533 if (updpar ==
nullptr) {
540 std::unique_ptr<const PseudoMeasurementOnTrack> newpseudo = std::make_unique<const PseudoMeasurementOnTrack>(
548 pseudostate->setMeasurement(std::move(newpseudo));
552 pseudostate->setMeasurementErrors(
errors);
553 pseudoupdated =
true;
564 *
track->perigeeParameters(),
575 if (
track !=
nullptr) {
576 track->info().addPatternReco(intrk1.
info());
577 track->info().addPatternReco(intrk2.
info());
588 const EventContext& ctx,
590 const Track & intrk1,
591 const Track & intrk2,
593 std::vector<MaterialEffectsOnTrack> & calomeots
595 ATH_MSG_DEBUG(
"--> entering GlobalChi2Fitter::mainCombinationStrategy");
600 const Track *indettrack = firstismuon ? &intrk2 : &intrk1;
601 const Track *muontrack = firstismuon ? &intrk1 : &intrk2;
603 auto [tmpfirstidpar, tmplastidpar] = getFirstLastIdPar(*indettrack);
604 std::unique_ptr<const TrackParameters> firstidpar =
unique_clone(tmpfirstidpar);
605 std::unique_ptr<const TrackParameters> lastidpar =
unique_clone(tmplastidpar);
607 if ((firstidpar ==
nullptr) || (lastidpar ==
nullptr)) {
621 std::unique_ptr<const TrackParameters> tp_closestmuon =
nullptr;
623 while (closestmuonmeas ==
nullptr) {
624 closestmuonmeas =
nullptr;
627 if ((**tsosit).measurementOnTrack() !=
nullptr) {
628 closestmuonmeas = (**tsosit).measurementOnTrack();
630 if (thispar !=
nullptr) {
631 const AmgVector(5) & parvec = thispar->parameters();
633 parvec[0], parvec[1], parvec[2], parvec[3], parvec[4], std::nullopt
647 std::unique_ptr<const TrackParameters> tmppar;
661 if ((tp_closestmuon !=
nullptr) && (cache.
m_msEntrance !=
nullptr)) {
666 std::unique_ptr<const std::vector<const TrackStateOnSurface *>> matvec;
668 if (tmppar !=
nullptr) {
669 const Surface & associatedSurface = tmppar->associatedSurface();
670 std::unique_ptr<Surface> muonsurf =
nullptr;
676 double radius = cylbounds->
r();
678 muonsurf = std::make_unique<CylinderSurface>(trans,
radius + 1, hlength);
683 associatedSurface.
center().z() > 0 ?
684 associatedSurface.
center().z() + 1 :
685 associatedSurface.
center().z() - 1
689 associatedSurface.
center().x(),
690 associatedSurface.
center().y(),
694 trans.translation() << newpos;
697 double rmin = discbounds->
rMin();
698 double rmax = discbounds->
rMax();
699 muonsurf = std::make_unique<DiscSurface>(trans, rmin, rmax);
703 if (muonsurf !=
nullptr) {
715 std::vector<const TrackStateOnSurface *> tmp_matvec;
717 if ((matvec !=
nullptr) && !matvec->empty()) {
718 tmp_matvec = *matvec;
719 delete tmp_matvec.back();
720 tmp_matvec.pop_back();
722 for (
auto &
i : tmp_matvec) {
741 if (tmppar ==
nullptr) {
755 if (tmppar ==
nullptr) {
759 AmgVector(5) newpars = tmppar->parameters();
765 double newp2 = oldp * oldp + (!firstismuon ? 2 : -2) * de * std::sqrt(
mass *
mass + oldp * oldp) + de * de;
772 tp_closestmuon=tmppar->associatedSurface().createUniqueTrackParameters(
773 newpars[0], newpars[1], newpars[2], newpars[3], newpars[4], std::nullopt
789 for (; itStates != endState; ++itStates) {
796 if ((*itStates)->materialEffectsOnTrack() !=
nullptr) {
803 const auto *
const pBaseMEOT = (*itStates)->materialEffectsOnTrack();
807 const auto *
const pMEOT =
static_cast<const MaterialEffectsOnTrack *
>((*itStates)->materialEffectsOnTrack());
808 if ((pMEOT->scatteringAngles() ==
nullptr) or (pMEOT->energyLoss() ==
nullptr)) {
822 trajectory.
trackStates().back()->setTrackParameters(
nullptr);
825 std::unique_ptr<const TrackParameters> firstscatpar;
826 std::unique_ptr<const TrackParameters> lastscatpar;
829 double newqoverpid = 0;
832 double de = std::abs(calomeots[1].energyLoss()->deltaE());
833 double sigmade = std::abs(calomeots[1].energyLoss()->sigmaDeltaE());
835 double pbefore = std::abs(1 / firstidpar->parameters()[
Trk::qOverP]);
836 double pafter = std::abs(1 / tp_closestmuon->parameters()[
Trk::qOverP]);
837 double elosspull = (pbefore - pafter - de) / sigmade;
839 if (std::abs(elosspull) > 10) {
840 if (elosspull > 10) {
841 newqoverpid = 1 / (de + pafter + 10 * sigmade);
843 newqoverpid = 1 / (de + pafter - 10 * sigmade);
846 if (tp_closestmuon->parameters()[
Trk::qOverP] * newqoverpid < 0) {
850 const AmgVector(5) & newpar = firstidpar->parameters();
851 firstidpar=firstidpar->associatedSurface().createUniqueTrackParameters(
852 newpar[0], newpar[1], newpar[2], newpar[3], newqoverpid, std::nullopt
860 if (lastidpar ==
nullptr) {
866 *(firstismuon ? tp_closestmuon.get() : lastidpar.get()),
867 calomeots[0].associatedSurface(),
874 if (firstscatpar ==
nullptr) {
880 *(firstismuon ? firstidpar : tp_closestmuon),
881 calomeots[2].associatedSurface(),
888 if (lastscatpar ==
nullptr) {
892 std::optional<TransportJacobian> jac1, jac2;
893 std::unique_ptr<const TrackParameters> elosspar;
895 double firstscatphi = 0;
896 double secondscatphi = 0;
897 double firstscattheta = 0;
898 double secondscattheta = 0;
899 double muonscatphi = 0;
900 double muonscattheta = 0;
902 const TrackParameters *idscatpar = !firstismuon ? firstscatpar.get() : lastscatpar.get();
903 const TrackParameters *muonscatpar = firstismuon ? firstscatpar.get() : lastscatpar.get();
905 newqoverpid = idscatpar->parameters()[
Trk::qOverP];
907 Amg::Vector3D calosegment = lastscatpar->position() - firstscatpar->position();
908 muonscatphi = calosegment.phi() - muonscatpar->parameters()[
Trk::phi];
910 if (std::abs(std::abs(muonscatphi) - 2 *
M_PI) < std::abs(muonscatphi)) {
911 muonscatphi += (muonscatphi < 0 ? 2 *
M_PI : -2 *
M_PI);
914 muonscattheta = calosegment.theta() - muonscatpar->parameters()[
Trk::theta];
917 for (
int i = 0;
i < 2;
i++) {
918 std::unique_ptr<const TrackParameters> tmpelosspar;
919 AmgVector(5) params1 = muonscatpar->parameters();
928 params1[0], params1[1], params1[2], params1[3], params1[4], std::nullopt
949 calomeots[1].associatedSurface(),
955 if ((tmpelosspar ==
nullptr) || (jac1 == std::nullopt)) {
964 const AmgVector(5) & newpars = tmpelosspar->parameters();
965 std::unique_ptr<const TrackParameters> elosspar2(tmpelosspar->associatedSurface().createUniqueTrackParameters(
966 newpars[0], newpars[1], newpars[2], newpars[3], newqoverpid, std::nullopt
970 elosspar = std::move(tmpelosspar);
973 std::unique_ptr<const TrackParameters> scat2(
m_propagator->propagateParameters(
977 calomeots[0].associatedSurface() :
978 calomeots[2].associatedSurface(),
991 calomeots[0].associatedSurface() :
992 calomeots[2].associatedSurface(),
1000 if ((scat2 ==
nullptr) || (jac2 == std::nullopt)) {
1005 for (
int j = 0; j < 5; j++) {
1006 for (
int k = 0;
k < 5;
k++) {
1008 for (
int l = 0;
l < 5;
l++) {
1009 jac3[j][
k] += (*jac2) (j,
l) * (*jac1) (
l,
k);
1018 jac4(0, 0) = jac3[0][2];
1019 jac4(1, 1) = jac3[1][3];
1020 jac4(0, 1) = jac3[0][3];
1021 jac4(1, 0) = jac3[1][2];
1023 jac4 = jac4.inverse();
1035 discsurf =
static_cast<const Trk::DiscSurface *
>(&scat2->associatedSurface());
1037 if (cylsurf !=
nullptr) {
1039 if (std::abs(std::abs(dloc1) -
length) < std::abs(dloc1)) {
1048 if (discsurf !=
nullptr) {
1049 if (std::abs(std::abs(dloc2) - 2 *
M_PI) < std::abs(dloc2)) {
1058 double dphi = jac4(0, 0) * dloc1 + jac4(0, 1) * dloc2;
1059 double dtheta = jac4(1, 0) * dloc1 + jac4(1, 1) * dloc2;
1065 muonscatphi += dphi;
1066 muonscattheta += dtheta;
1068 double idscatphi = idscatpar->parameters()[
Trk::phi] - (scat2->parameters()[
Trk::phi] + dphi);
1070 if (std::abs(std::abs(idscatphi) - 2 *
M_PI) < std::abs(idscatphi)) {
1071 idscatphi += ((idscatphi < 0) ? 2 *
M_PI : -2 *
M_PI);
1074 double idscattheta = idscatpar->parameters()[
Trk::theta] - (scat2->parameters()[
Trk::theta] + dtheta);
1077 firstscatphi = muonscatphi;
1078 secondscatphi = idscatphi;
1079 firstscattheta = muonscattheta;
1080 secondscattheta = idscattheta;
1082 firstscatphi = -idscatphi;
1083 secondscatphi = -muonscatphi;
1084 firstscattheta = -idscattheta;
1085 secondscattheta = -muonscattheta;
1089 AmgVector(5) params2 = scat2->parameters();
1097 firstscatpar=scat2->associatedSurface().createUniqueTrackParameters(
1098 params2[0], params2[1], params2[2], params2[3], params2[4], std::nullopt
1100 idscatpar = firstscatpar.get();
1108 if (startPar !=
nullptr) {
1117 rot.col(2) = trackdir;
1120 trans.linear().matrix() << rot;
1121 trans.translation() << startPar->position() - .1 * trackdir;
1132 if (curvlinpar !=
nullptr) {
1133 startPar.reset(curvlinpar);
1137 firstscatpar = std::move(scat2);
1141 std::unique_ptr<GXFMaterialEffects> firstscatmeff = std::make_unique<GXFMaterialEffects>(calomeots[0]);
1142 std::unique_ptr<GXFMaterialEffects> elossmeff = std::make_unique<GXFMaterialEffects>(calomeots[1]);
1143 std::unique_ptr<GXFMaterialEffects> secondscatmeff = std::make_unique<GXFMaterialEffects>(calomeots[2]);
1145 double pull1 = std::abs(firstscatphi / firstscatmeff->sigmaDeltaPhi());
1146 double pull2 = std::abs(secondscatphi / secondscatmeff->sigmaDeltaPhi());
1149 for (
auto &
i : tmp_matvec) {
1154 firstscatmeff->setScatteringAngles(firstscatphi, firstscattheta);
1155 secondscatmeff->setScatteringAngles(secondscatphi, secondscattheta);
1158 elossmeff->setdelta_p(1000 * (lastscatpar->parameters()[
Trk::qOverP] - newqoverpid));
1160 elossmeff->setdelta_p(1000 * (newqoverpid - firstscatpar->parameters()[
Trk::qOverP]));
1163 elossmeff->setSigmaDeltaE(calomeots[1].energyLoss()->sigmaDeltaE());
1165 trajectory.
addMaterialState(std::make_unique<GXFTrackState>(std::move(firstscatmeff), std::move(firstscatpar)), -1);
1166 trajectory.
addMaterialState(std::make_unique<GXFTrackState>(std::move(elossmeff), std::move(elosspar)), -1);
1167 trajectory.
addMaterialState(std::make_unique<GXFTrackState>(std::move(secondscatmeff), std::move(lastscatpar)), -1);
1170 for (
auto &
i : tmp_matvec) {
1177 if (startPar ==
nullptr) {
1182 (pull1 > 5 || pull2 > 5) &&
1188 bool largegap =
false;
1189 double previousz = 0;
1191 for (itStates2 = beginStates2; itStates2 != endState2; ++itStates2) {
1196 if (meff !=
nullptr) {
1202 if ( mefot and mefot->energyLoss() and
1203 std::abs(mefot->energyLoss()->deltaE()) > 250 &&
1204 mefot->energyLoss()->sigmaDeltaE() < 1.e-9
1217 !(itStates2 == beginStates2 || itStates2 == beginStates2 + 1) &&
1226 itStates2 == endState2 - 1 &&
1229 std::abs(tpar->
position().z()) < 13000
1231 std::unique_ptr<const TrackParameters> pseudopar(tpar->
clone());
1235 std::unique_ptr<const PseudoMeasurementOnTrack> newpseudo = std::make_unique<const PseudoMeasurementOnTrack>(
1238 pseudopar->associatedSurface()
1241 std::unique_ptr<GXFTrackState> pseudostate = std::make_unique<GXFTrackState>(std::move(newpseudo), std::move(pseudopar));
1248 pseudostate->setMeasurementErrors(
errors);
1255 std::abs(trajectory.
trackStates().back()->position().z()) > 20000 &&
1256 std::abs(previousz) < 12000
1262 previousz = trajectory.
trackStates().back()->position().z();
1280 const EventContext& ctx,
1282 const Track & intrk1,
1283 const Track & intrk2,
1285 std::vector<MaterialEffectsOnTrack> & calomeots
1287 ATH_MSG_DEBUG(
"--> entering GlobalChi2Fitter::backupCombinationStrategy");
1289 bool firstismuon =
false;
1290 const Track *indettrack = &intrk1;
1293 indettrack = &intrk2;
1307 if (pParametersVector->size() > 1)
1308 firstidpar = (*pParametersVector)[1];
1310 firstidpar = pParametersVector->back();
1312 std::unique_ptr<const TrackParameters> lastidpar =
nullptr;
1313 if ((firstidpar !=
nullptr) && (cache.
m_caloEntrance !=
nullptr))
1317 if (lastidpar ==
nullptr) {
1318 lastidpar.reset(pParametersVector->back()->clone());
1321 std::unique_ptr < const TrackParameters > firstscatpar;
1322 std::unique_ptr < const TrackParameters > lastscatpar;
1323 std::unique_ptr < const TrackParameters > elosspar;
1328 lastidpar->position(),
1329 lastidpar->momentum(),
1331 lastidpar->position()
1338 calomeots[0].associatedSurface(),
1344 if (!firstscatpar) {
1348 std::unique_ptr<const TrackParameters> tmppar(
1352 calomeots[1].associatedSurface(),
1365 double oldp = std::abs(1 / tmppar->parameters()[
Trk::qOverP]);
1366 double de = std::abs(calomeots[1].energyLoss()->deltaE());
1368 double newp2 = oldp * oldp - 2 * de * std::sqrt(
mass *
mass + oldp * oldp) + de * de;
1374 double newqoverp =
sign / std::sqrt(newp2);
1379 tmppar->associatedSurface().createUniqueTrackParameters(
1386 calomeots[2].associatedSurface(),
1399 calomeots[2].associatedSurface(),
1413 calomeots[1].associatedSurface(),
1425 double newqoverp =
sign /
1426 (1. / std::abs(elosspar->parameters()[
Trk::qOverP]) +
1427 std::abs(calomeots[1].energyLoss()->deltaE()));
1431 std::unique_ptr<const TrackParameters>tmppar(
1432 elosspar->associatedSurface().createUniqueTrackParameters(
1440 calomeots[0].associatedSurface(),
1447 if (!firstscatpar) {
1452 for (; itStates != endState; ++itStates) {
1460 if ((*itStates)->materialEffectsOnTrack() !=
nullptr) {
1473 std::unique_ptr<GXFMaterialEffects> firstscatmeff = std::make_unique<GXFMaterialEffects>(calomeots[0]);
1474 std::unique_ptr<GXFMaterialEffects> elossmeff = std::make_unique<GXFMaterialEffects>(calomeots[1]);
1475 std::unique_ptr<GXFMaterialEffects> secondscatmeff = std::make_unique<GXFMaterialEffects>(calomeots[2]);
1479 sigmadp = calomeots[1].energyLoss()->sigmaDeltaE();
1480 elossmeff->setSigmaDeltaE(sigmadp);
1483 elossmeff->setdelta_p(
dp);
1485 trajectory.
addMaterialState(std::make_unique<GXFTrackState>(std::move(firstscatmeff), std::move(firstscatpar)), -1);
1486 trajectory.
addMaterialState(std::make_unique<GXFTrackState>(std::move(elossmeff), std::move(elosspar)), -1);
1487 trajectory.
addMaterialState(std::make_unique<GXFTrackState>(std::move(secondscatmeff), std::move(lastscatpar)), -1);
1490 const Surface *triggersurf1 =
nullptr;
1491 const Surface *triggersurf2 =
nullptr;
1495 bool seenmdt =
false;
1496 bool mdtbetweenphihits =
false;
1500 itStates2 = (!firstismuon ? beginStates2 : endState - 1);
1501 itStates2 != (!firstismuon ? endState2 : beginStates - 1);
1502 (!firstismuon ? ++itStates2 : --itStates2)
1505 ((*itStates2)->measurementOnTrack() ==
nullptr) ||
1510 const auto *
const pMeasurement = (*itStates2)->measurementOnTrack();
1511 const Surface *surf = &pMeasurement->associatedSurface();
1515 if (isCompetingRIOsOnTrack) {
1517 rot = &crot->rioOnTrack(0);
1520 rot =
static_cast<const RIO_OnTrack *
>(pMeasurement);
1527 (rot !=
nullptr) && (
1533 bool measphi =
true;
1537 if (std::abs(dotprod1) > .5 || std::abs(dotprod2) > .5) {
1543 (*itStates2)->trackParameters() !=
nullptr ?
1544 (*itStates2)->trackParameters()->position() :
1546 if (triggersurf1 !=
nullptr) {
1547 triggerpos2 = thispos;
1548 triggersurf2 = surf;
1550 mdtbetweenphihits =
true;
1553 triggerpos1 = thispos;
1554 triggersurf1 = surf;
1560 double mdttrig1 = 999999;
1561 double mdttrig2 = 999999;
1562 const Surface *mdtsurf1 =
nullptr;
1563 const Surface *mdtsurf2 =
nullptr;
1566 itStates2 = (!firstismuon ? beginStates2 : endState - 1);
1567 itStates2 != (!firstismuon ? endState2 : beginStates - 1);
1568 (!firstismuon ? ++itStates2 : --itStates2)
1570 const Surface *surf =
nullptr;
1572 ((*itStates2)->measurementOnTrack() !=
nullptr) &&
1575 surf = &(*itStates2)->measurementOnTrack()->associatedSurface();
1578 if (surf ==
nullptr) {
1581 const auto *
const pThisMeasurement = (*itStates2)->measurementOnTrack();
1585 if (isCompetingRioOnTrack) {
1587 rot = &crot->rioOnTrack(0);
1590 rot =
static_cast<const RIO_OnTrack *
>(pThisMeasurement);
1596 (*itStates2)->trackParameters() !=
nullptr ?
1597 (*itStates2)->trackParameters()->position() :
1598 pThisMeasurement->globalPosition();
1599 if (triggerpos1.mag() > 1 && (globpos - triggerpos1).
mag() < mdttrig1) {
1600 mdttrig1 = (globpos - triggerpos1).
mag();
1603 if (triggerpos2.mag() > 1 && (globpos - triggerpos2).
mag() < mdttrig2) {
1604 mdttrig2 = (globpos - triggerpos2).
mag();
1611 std::vector<GXFTrackState *> outlierstates;
1612 std::vector<GXFTrackState *> outlierstates2;
1614 outlierstates.reserve(10);
1615 outlierstates2.reserve(10);
1617 std::unique_ptr<PseudoMeasurementOnTrack> newpseudo;
1619 for (itStates2 = beginStates2; itStates2 != endState2; ++itStates2) {
1620 const auto *
const pMeasurement{(*itStates2)->measurementOnTrack()};
1622 bool isStraightLine =
1623 pMeasurement !=
nullptr ?
1630 (newpseudo ==
nullptr) && (
1631 (itStates2 == beginStates2 || itStates2 == beginStates2 + 1) &&
1632 std::abs(pMeasurement->globalPosition().z()) < 10000
1635 std::unique_ptr<const TrackParameters> par2;
1636 if (((*itStates2)->trackParameters() !=
nullptr) && nphi > 99) {
1637 par2.reset((*itStates2)->trackParameters()->clone());
1649 if (par2 ==
nullptr) {
1655 newpseudo = std::make_unique<PseudoMeasurementOnTrack>(
1658 par2->associatedSurface()
1661 std::unique_ptr<GXFTrackState> firstpseudo = std::make_unique<GXFTrackState>(std::move(newpseudo), std::move(par2));
1668 firstpseudo->setMeasurementErrors(
errors);
1669 firstpseudostate = firstpseudo.get();
1675 if (isPseudo && !firstismuon) {
1679 if ((**itStates2).materialEffectsOnTrack() !=
nullptr) {
1689 ((**itStates2).measurementOnTrack() !=
nullptr) &&
1690 &(**itStates2).measurementOnTrack()->associatedSurface() == triggersurf1 &&
1691 (mdtsurf1 !=
nullptr)
1693 std::unique_ptr<Amg::Transform3D> transf = std::make_unique<Amg::Transform3D>(mdtsurf1->
transform());
1695 transf->translation() << triggerpos1;
1700 std::unique_ptr<const PseudoMeasurementOnTrack> newpseudo = std::make_unique<const PseudoMeasurementOnTrack>(
1704 std::unique_ptr<GXFTrackState> pseudostate1 = std::make_unique<GXFTrackState>(std::move(newpseudo),
nullptr);
1711 pseudostate1->setMeasurementErrors(
errors);
1712 outlierstates2.push_back(pseudostate1.get());
1717 ((**itStates2).measurementOnTrack() !=
nullptr) &&
1718 &(**itStates2).measurementOnTrack()->associatedSurface() == triggersurf2 &&
1719 mdtbetweenphihits &&
1720 (mdtsurf2 !=
nullptr)
1722 std::unique_ptr<Amg::Transform3D> transf = std::make_unique<Amg::Transform3D>(mdtsurf2->
transform());
1723 transf->translation() << triggerpos2;
1728 std::unique_ptr<const PseudoMeasurementOnTrack> newpseudo = std::make_unique<const PseudoMeasurementOnTrack>(
1732 std::unique_ptr<GXFTrackState> pseudostate2 = std::make_unique<GXFTrackState>(std::move(newpseudo),
nullptr);
1739 pseudostate2->setMeasurementErrors(
errors);
1741 outlierstates2.push_back(pseudostate2.get());
1757 outlierstates.push_back(trajectory.
trackStates().back().get());
1773 std::unique_ptr<Trk::Track> tmp_track(
1774 myfit(ctx, cache, trajectory, *startpar2,
false,
muon));
1781 std::abs(trajectory.
residuals().tail<1>()(0) / trajectory.
errors().tail<1>()(0)) > 10
1787 if (firstpseudostate !=
nullptr) {
1792 std::unique_ptr<const PseudoMeasurementOnTrack> newpseudo = std::make_unique<const PseudoMeasurementOnTrack>(
1801 for (
int j = 0; j < (
int) trajectory.
trackStates().size(); j++) {
1802 for (
const auto &
i : outlierstates2) {
1808 for (
const auto &
i : outlierstates) {
1816 itStates = (firstismuon ? beginStates2 : endState - 1);
1817 itStates != (firstismuon ? endState2 : beginStates - 1);
1818 (firstismuon ? ++itStates : --itStates)
1824 makeProtoState(cache, trajectory, *itStates, (firstismuon ? -1 : 0));
1837 std::unique_ptr<Track>
1839 const EventContext& ctx,
1840 const Track& inputTrack,
1844 ATH_MSG_DEBUG(
"--> entering GlobalChi2Fitter::fit(Track,)");
1857 return std::unique_ptr<Track>(
1858 fitIm(ctx, cache, inputTrack, runOutlier, matEffects));
1863 const Track &inputTrack,
1868 const EventContext& ctx = Gaudi::Hive::currentContext();
1872 delete alignCache.m_derivMatrix;
1873 alignCache.m_derivMatrix =
nullptr;
1875 delete alignCache.m_fullCovarianceMatrix;
1876 alignCache.m_fullCovarianceMatrix =
nullptr;
1877 alignCache.m_iterationsOfLastFit = 0;
1880 fitIm(ctx, cache, inputTrack, runOutlier, matEffects);
1881 if(newTrack !=
nullptr){
1886 alignCache.m_iterationsOfLastFit = cache.
m_lastiter;
1893 const EventContext& ctx,
1895 const Track& inputTrack,
1900 ATH_MSG_DEBUG(
"--> entering GlobalChi2Fitter::fit(Track,,)");
1916 ATH_MSG_WARNING(
"Track without track parameters, cannot perform fit");
1924 if (minpar ==
nullptr) {
1942 const Surface *firsthitsurf =
nullptr;
1943 const Surface *lasthitsurf =
nullptr;
1945 bool hasmuon =
false;
1946 bool iscombined =
false;
1947 bool seenphimeas =
false;
1951 for (; itStates != endState; ++itStates) {
1952 const auto *
const pMeasurement = (**itStates).measurementOnTrack();
1954 (pMeasurement ==
nullptr) &&
1955 ((**itStates).materialEffectsOnTrack() !=
nullptr) &&
1961 if (pMeasurement !=
nullptr) {
1962 const Surface *surf = &pMeasurement->associatedSurface();
1971 if (firsthitsurf ==
nullptr) {
1972 firsthitsurf = surf;
1980 if ((**itStates).trackParameters() !=
nullptr) {
1981 lastidpar = (**itStates).trackParameters();
1982 if (firstidpar ==
nullptr) {
1983 firstidpar = lastidpar;
1991 if (std::abs(dotprod1) < .5 && std::abs(dotprod2) < .5) {
1993 if (std::abs(surf->
center().z()) > 13000) {
1996 if (surf->
center().perp() > 9000 && std::abs(surf->
center().z()) < 13000) {
2008 if (iscombined && seenphimeas && (phiem || phibo)) {
2024 if (firstidpar == lastidpar) {
2025 firstidpar = lastidpar =
nullptr;
2035 (firstidpar !=
nullptr)
2047 std::unique_ptr<Track> tmptrack =
nullptr;
2058 tmptrack.reset(
myfit(ctx, cache, trajectory, *minpar,
false, matEffects));
2061 if (tmptrack ==
nullptr) {
2070 bool isbrem =
false;
2072 unsigned int n_brem=0;
2074 for (std::unique_ptr<GXFTrackState> & state : trajectory.
trackStates()) {
2077 if (meff !=
nullptr) {
2085 std::optional<Amg::Vector2D> locpos(state->associatedSurface().globalToLocal(layerpars->
position()));
2086 const Amg::Vector3D layerNormal(state->associatedSurface().normal(*locpos));
2087 double costracksurf = 1.;
2089 costracksurf = std::abs(layerNormal.dot(layerpars->
momentum().unit()));
2091 double oldde = meff->
deltaE();
2093 std::unique_ptr<EnergyLoss> eloss;
2094 double sigmascat = 0;
2096 if (matprop !=
nullptr) {
2097 eloss = std::make_unique<EnergyLoss>(
2098 m_elosstool->energyLoss(*matprop,
p, 1. / costracksurf,
2100 sigmascat = std::sqrt(
m_scattool->sigmaSquare(*matprop,
p, 1. / costracksurf, matEffects));
2102 if (eloss !=
nullptr) {
2107 sigmascat = std::sqrt(
m_scattool->sigmaSquare(tmpprop,
p, 1. / costracksurf, matEffects));
2125 }
else if (eloss !=
nullptr) {
2136 refpars[0], refpars[1], refpars[2], refpars[3], refpars[4], std::nullopt
2144 trajectory.
brems().clear();
2146 trajectory.
brems().resize(1);
2147 trajectory.
brems()[0] = bremdp;
2158 std::unique_ptr<Track>
track(
myfit(ctx, cache, trajectory, *minpar, runOutlier, matEffects));
2160 bool pseudoupdated =
false;
2162 if ((
track !=
nullptr) && hasid && hasmuon) {
2163 for (std::unique_ptr<GXFTrackState> & pseudostate : trajectory.
trackStates()) {
2165 (pseudostate ==
nullptr) ||
2167 pseudostate->fitQuality().chiSquared() < 10
2173 std::unique_ptr<const TrackParameters> updpar(
m_updator->removeFromState(
2175 pseudostate->measurement()->localParameters(),
2176 pseudostate->measurement()->localCovariance()
2179 if (updpar ==
nullptr) {
2186 std::unique_ptr<const PseudoMeasurementOnTrack> newpseudo = std::make_unique<const PseudoMeasurementOnTrack>(
2192 pseudostate->setMeasurement(std::move(newpseudo));
2196 pseudostate->setMeasurementErrors(
errors);
2197 pseudoupdated =
true;
2200 if (pseudoupdated) {
2213 if (
track !=
nullptr) {
2216 track->info().addPatternReco(old_info);
2219 return track.release();
2222 std::unique_ptr<Track>
2224 const EventContext& ctx,
2230 ATH_MSG_DEBUG(
"--> entering GlobalChi2Fitter::fit(PRDS,TP,)");
2233 for (
const auto *prd : prds) {
2234 const Surface & prdsurf = (*prd).detectorElement()->surface((*prd).identify());
2239 plsurf =
static_cast < const PlaneSurface *
>(&prdsurf);
2246 if ((slsurf ==
nullptr) && (plsurf ==
nullptr)) {
2247 ATH_MSG_ERROR(
"Surface is neither PlaneSurface nor StraightLineSurface!");
2252 }
else if (slsurf !=
nullptr) {
2261 }
else if (plsurf !=
nullptr) {
2262 if (param.covariance() !=
nullptr) {
2284 if (rot !=
nullptr) {
2285 rots.push_back(rot);
2289 std::unique_ptr<Track>
track =
2290 fit(ctx, rots, param, runOutlier, matEffects);
2292 for (
const auto *rot : rots) {
2299 std::unique_ptr<Track>
2301 const EventContext& ctx,
2302 const Track& inputTrack,
2307 ATH_MSG_DEBUG(
"--> entering GlobalChi2Fitter::fit(Track,Meas'BaseSet,,)");
2322 if (minpar ==
nullptr) {
2340 for (; itStates != endState; ++itStates) {
2345 (trajectory.
trackStates().back()->materialEffects() !=
nullptr) &&
2346 trajectory.
trackStates().back()->materialEffects()->sigmaDeltaE() > 50.001
2348 trajectory.
trackStates().back()->materialEffects()->setKink(
true);
2353 MeasurementSet::const_iterator itSet = addMeasColl.begin();
2354 MeasurementSet::const_iterator itSetEnd = addMeasColl.end();
2356 for (; itSet != itSetEnd; ++itSet) {
2357 if ((*itSet) ==
nullptr) {
2358 ATH_MSG_WARNING(
"There is an empty MeasurementBase object in the track! Skip this object..");
2365 std::unique_ptr<Track>
track(
myfit(ctx, cache, trajectory, *minpar, runOutlier, matEffects));
2368 if (
track !=
nullptr) {
2375 track->fitQuality()->numberDoF() != 0 ?
2376 track->fitQuality()->chiSquared() /
track->fitQuality()->numberDoF() :
2379 if (
m_extensioncuts && runOutlier && newqual > 2 && newqual > 2 * oldqual) {
2383 track.reset(
nullptr);
2387 if (
track !=
nullptr) {
2396 std::unique_ptr<Track>
2398 const EventContext& ctx,
2404 ATH_MSG_DEBUG(
"--> entering GlobalChi2Fitter::fit(Track,PRDS,)");
2408 for (
const auto *prd : prds) {
2409 const Surface & prdsurf = (*prd).detectorElement()->surface((*prd).identify());
2412 std::unique_ptr<const TrackParameters>trackparForCorrect(
2428 rot =
m_ROTcreator->correct(*prd, *trackparForCorrect, ctx);
2431 if (rot !=
nullptr) {
2432 rots.push_back(rot);
2436 std::unique_ptr<Track>
track =
fit(ctx,intrk, rots, runOutlier, matEffects);
2438 for (
const auto *rot : rots) {
2446 const EventContext& ctx,
2452 ATH_MSG_DEBUG(
"--> entering GlobalChi2Fitter::fit(Meas'BaseSet,,)");
2465 for (
const auto *itSet : rots) {
2466 if (itSet ==
nullptr) {
2467 ATH_MSG_WARNING(
"There is an empty MeasurementBase object in the track! Skip this object..");
2473 std::unique_ptr<const TrackParameters> startpar(param.
clone());
2476 matEffects ==
muon &&
2482 myfit(ctx,cache, trajectory, *startpar, runOutlier, matEffects);
2500 std::unique_ptr<const PseudoMeasurementOnTrack> newpseudo = std::make_unique<const PseudoMeasurementOnTrack>(
2506 trajectory.
trackStates().front()->setMeasurement(std::move(newpseudo));
2513 std::unique_ptr<const PseudoMeasurementOnTrack> newpseudo = std::make_unique<const PseudoMeasurementOnTrack>(
2519 trajectory.
trackStates().back()->setMeasurement(std::move(newpseudo));
2526 refpars[0], refpars[1], refpars[2], refpars[3], refpars[4], std::nullopt
2531 myfit(ctx,cache, trajectory, *startpar, runOutlier, matEffects);
2542 refpars[0], refpars[1], refpars[2], refpars[3], refpars[4], std::nullopt
2549 firstpar = trajectory.
trackStates().front()->trackParameters();
2554 std::unique_ptr<const PseudoMeasurementOnTrack> newpseudo = std::make_unique<const PseudoMeasurementOnTrack>(
2560 trajectory.
trackStates().front()->setMeasurement(std::move(newpseudo));
2568 lastpar = trajectory.
trackStates().back()->trackParameters();
2573 std::unique_ptr<const PseudoMeasurementOnTrack> newpseudo = std::make_unique<const PseudoMeasurementOnTrack>(
2579 trajectory.
trackStates().back()->setMeasurement(std::move(newpseudo));
2587 std::unique_ptr<Track>
track;
2589 if (startpar !=
nullptr) {
2590 track.reset(
myfit(ctx,cache, trajectory, *startpar, runOutlier, matEffects));
2593 if (
track !=
nullptr) {
2623 std::unique_ptr<GXFMaterialEffects> newmeff;
2631 newmeff = std::make_unique<GXFMaterialEffects>(*meff);
2635 double sigmascat = std::sqrt(
m_scattool->sigmaSquare(
2656 newmeff = std::make_unique<GXFMaterialEffects>(newmeot);
2672 ((**trajectory.
trackStates().rbegin()).trackParameters() !=
nullptr)
2674 double delta_p = 1000 * (
2676 (**trajectory.
trackStates().rbegin()).trackParameters()->
2680 newmeff->setdelta_p(delta_p);
2691 bool isoutlier =
false;
2720 seg =
static_cast<const Segment *
>(measbase);
2730 for (
int i = 0;
i <
imax;
i++) {
2733 std::unique_ptr<GXFTrackState> ptsos = std::make_unique<GXFTrackState>(
2734 std::unique_ptr<const MeasurementBase>(measbase2->
clone()),
2735 std::unique_ptr<const TrackParameters>(newtrackpar !=
nullptr ? newtrackpar->
clone() :
nullptr)
2738 double sinstereo = 0;
2751 bool measphi =
false;
2754 bool rotated =
false;
2792 const double diagonalProduct =
covmat(0, 0) *
covmat(1, 1);
2793 const double element01Sq =
covmat(0, 1) *
covmat(0, 1);
2794 const double sqrtTerm = std::sqrt(
2795 (traceCov) * (traceCov) - 4. * (diagonalProduct - element01Sq)
2801 sinstereo =
std::sin(0.5 * std::asin(2 *
covmat(0, 1) / (-sqrtTerm)));
2819 if (std::abs(dotprod1) < .5 && std::abs(dotprod2) < .5) {
2827 int param_index = 0;
2834 errors[1] = std::sqrt(
covmat(param_index, param_index));
2839 errors[2] = std::sqrt(
covmat(param_index, param_index));
2844 errors[3] = std::sqrt(
covmat(param_index, param_index));
2849 errors[4] = std::sqrt(
covmat(param_index, param_index));
2870 ptsos->setMeasurementErrors(
errors);
2871 ptsos->setSinStereo(sinstereo);
2872 ptsos->setMeasurementType(hittype);
2873 ptsos->setMeasuresPhi(measphi);
2894 if (tvol ==
nullptr) {
2901 if (confinedLayers !=
nullptr) {
2906 for (; layerIter != layerVector.end(); ++layerIter) {
2908 if (*layerIter !=
nullptr) {
2913 if ((layIndex.
value() == 0) || ((*layerIter)->layerMaterialProperties() ==
nullptr)) {
2924 disclay =
static_cast<const DiscLayer *
>((*layerIter));
2926 if (disclay !=
nullptr) {
2927 if (disclay->
center().z() < 0) {
2932 }
else if (cyllay !=
nullptr) {
2943 for (
size_t ib = 0 ;
ib < bsurf.size(); ++
ib) {
2944 const Layer *
layer = bsurf[
ib]->surfaceRepresentation().materialLayer();
2946 if (
layer ==
nullptr)
continue;
2950 if ((layIndex.
value() == 0) || (
layer->layerMaterialProperties() ==
nullptr)) {
2962 discsurf =
static_cast<const DiscSurface *
>(&
layer->surfaceRepresentation());
2964 if (discsurf !=
nullptr) {
2966 discsurf->
center().z() < 0 &&
2971 discsurf->
center().z() > 0 &&
2977 (cylsurf !=
nullptr) &&
2983 if ((cylsurf ==
nullptr) && (discsurf ==
nullptr)) {
2990 if (confinedVolumes !=
nullptr) {
2996 for (; volIter != volIterEnd; ++volIter) {
2997 if (*volIter !=
nullptr) {
3013 std::optional<std::pair<Amg::Vector3D, double>>
3026 const double *
pos = parforextrap.
position().data();
3036 double r = (std::abs(currentqoverp) > 1
e-10) ? -sintheta / (currentqoverp * 300. *
field[2]) : 1e6;
3037 double xc = parforextrap.
position().x() -
r * sinphi;
3038 double yc = parforextrap.
position().y() +
r * cosphi;
3041 double delta_s = (surf.
center().z() -
z0) / costheta;
3053 double costracksurf = std::abs(costheta);
3055 return std::make_pair(
intersect, costracksurf);
3058 std::optional<std::pair<Amg::Vector3D, double>>
3073 const double *
pos = parforextrap.
position().data();
3081 double r = (std::abs(currentqoverp) > 1
e-10) ? -sintheta / (currentqoverp * 300. *
field[2]) : 1e6;
3082 double xc = parforextrap.
position().x() -
r * sinphi;
3083 double yc = parforextrap.
position().y() +
r * cosphi;
3086 double d = xc * xc + yc * yc;
3087 double rcyl = surf.
bounds().
r();
3088 double mysqrt = ((
r + rcyl) * (
r + rcyl) -
d) * (
d - (
r - rcyl) * (
r - rcyl));
3094 mysqrt = std::sqrt(mysqrt);
3095 double firstterm = xc / 2 + (xc * (rcyl * rcyl -
r *
r)) / (2 *
d);
3096 double secondterm = (mysqrt * yc) / (2 *
d);
3097 double x1 = firstterm + secondterm;
3098 double x2 = firstterm - secondterm;
3099 firstterm = yc / 2 + (yc * (rcyl * rcyl -
r *
r)) / (2 *
d);
3100 secondterm = (mysqrt * xc) / (2 *
d);
3101 double y1 = firstterm - secondterm;
3102 double y2 = firstterm + secondterm;
3105 double dist1 = (
x -
x1) * (
x -
x1) + (
y -
y1) * (
y -
y1);
3106 double dist2 = (
x -
x2) * (
x -
x2) + (
y -
y2) * (
y -
y2);
3108 if (dist1 < dist2) {
3116 double phi1 = std::atan2(
y - yc,
x - xc);
3117 double deltaphi = phi1 -
phi0;
3119 coercePhiCoordinateRange(deltaphi);
3121 double delta_z =
r * deltaphi / tantheta;
3122 double z =
z0 + delta_z;
3131 double phidir = parforextrap.parameters()[
Trk::phi] + deltaphi;
3132 coercePhiCoordinateRange(phidir);
3136 double costracksurf = std::abs(normal.unit().dot(trackdir));
3138 return std::make_pair(
intersect, costracksurf);
3145 std::vector<std::pair<const Layer *, const Layer *>> &
layers,
3150 std::vector<std::unique_ptr<GXFTrackState>> &
states = trajectory.
trackStates();
3151 std::vector<std::unique_ptr<GXFTrackState>> oldstates = std::move(
states);
3162 for (
int i = 0;
i <= indexoffset;
i++) {
3172 for (
int i = indexoffset + 1;
i < (
int) oldstates.size();
i++) {
3173 double rmeas = oldstates[
i]->position().perp();
3174 double zmeas = oldstates[
i]->position().z();
3182 while (layerindex < (
int)
layers.size()) {
3184 double costracksurf = 0.0;
3205 if (oldstates[
i]->trackParameters() !=
nullptr) {
3206 double rlayer = cylsurf->
bounds().
r();
3207 if (std::abs(rmeas - rlayer) < std::abs(parforextrap->
position().perp() - rlayer)) {
3208 parforextrap = oldstates[
i]->trackParameters();
3224 if (cylsurf->
bounds().
r() > rmeas)
break;
3234 if (oldstates[
i]->trackParameters() !=
nullptr) {
3235 double zlayer = discsurf->
center().z();
3236 if (std::abs(zmeas - zlayer) < std::abs(parforextrap->
position().z() - zlayer)) {
3237 parforextrap = oldstates[
i]->trackParameters();
3248 if (std::abs(discsurf->
center().z()) > std::abs(zmeas))
break;
3250 throw std::logic_error(
"Unhandled surface.");
3258 if (matprop ==
nullptr) {
3271 double actualx0 =
X0 / costracksurf;
3272 double de = -std::abs(
3276 (
m_p != 0.0 ? std::abs(
m_p) : std::abs(1. / currentqoverp)),
3279 double sigmascat = std::sqrt(
m_scattool->sigmaSquare(
3281 (
m_p != 0.0 ? std::abs(
m_p) : std::abs(1. / currentqoverp)),
3285 std::unique_ptr<GXFMaterialEffects> meff = std::make_unique<GXFMaterialEffects>();
3286 meff->setDeltaE(de);
3287 meff->setScatteringSigmas(sigmascat / sintheta, sigmascat);
3288 meff->setX0(actualx0);
3289 meff->setSurface(&
layer->surfaceRepresentation());
3290 meff->setMaterialProperties(matprop);
3296 std::unique_ptr<EnergyLoss> eloss;
3299 eloss = std::make_unique<EnergyLoss>(
m_elosstool->energyLoss(
3301 (
m_p != 0.0 ? std::abs(
m_p) : std::abs(1. / currentqoverp)),
3306 if (eloss !=
nullptr) {
3307 meff->setSigmaDeltaE(eloss->sigmaDeltaE());
3312 meff->setDeltaE(-5);
3314 meff->setScatteringSigmas(0, 0);
3317 meff->setSigmaDeltaE(50);
3318 if (eloss !=
nullptr) {
3319 meff->setSigmaDeltaEPos(eloss->sigmaPlusDeltaE());
3320 meff->setSigmaDeltaENeg(eloss->sigmaMinusDeltaE());
3325 "X0: " << meff->x0() <<
" qoverp: " << currentqoverp <<
3326 " sigmascat " << meff->sigmaDeltaTheta() <<
" eloss: " << meff->deltaE() <<
3327 " sigma eloss: " << meff->sigmaDeltaE()
3334 std::unique_ptr<GXFTrackState> matstate = std::make_unique<GXFTrackState>(
3336 std::unique_ptr<const TrackParameters>()
3354 std::vector<std::pair<const Layer *, const Layer *>> &
layers,
3355 std::vector<std::pair<const Layer *, const Layer *>> & upstreamlayers,
3356 const std::vector<std::unique_ptr<GXFTrackState>> & oldstates,
3365 upstreamlayers.reserve(5);
3384 double lastz = laststate->
position().z();
3385 double lastr = laststate->
position().perp();
3392 double slope = (tantheta != 0) ? 1 / tantheta : 0;
3398 std::vector < const Layer *>::const_iterator
it;
3399 std::vector < const Layer *>::const_iterator itend;
3417 for (;
it != itend; ++
it) {
3422 if (std::abs((*it)->surfaceRepresentation().center().z()) > std::abs(lastz)) {
3430 const DiscBounds *discbounds = (
const DiscBounds *) (&(*it)->surfaceRepresentation().bounds());
3435 if (discbounds->
rMax() < firstr || discbounds->
rMin() > lastr) {
3439 double rintersect = firstr + ((*it)->surfaceRepresentation().center().z() - firstz) / slope;
3442 rintersect < discbounds->rMin() - 50 ||
3443 rintersect > discbounds->
rMax() + 50
3453 if ((*
it) == endlayer) {
3465 std::abs((*it)->surfaceRepresentation().center().z()) < std::abs(firstz) ||
3468 upstreamlayers.emplace_back((
Layer *)
nullptr, (*it));
3475 (*
it) != startlayer &&
3476 (std::abs((*it)->surfaceRepresentation().center().z()) > std::abs(firstz2) ||
3477 (*it) == startlayer2)
3491 if (barrelcylinder->surfaceRepresentation().bounds().r() > lastr) {
3498 double zintersect = firstz + (barrelcylinder->surfaceRepresentation().bounds().r() - firstr) * slope;
3500 if (std::abs(zintersect - barrelcylinder->surfaceRepresentation().center().z()) >
3501 ((
const CylinderSurface*)(&barrelcylinder->surfaceRepresentation()))->bounds().halflengthZ() + 50) {
3505 if (barrelcylinder == endlayer) {
3512 if (barrelcylinder->surfaceRepresentation().bounds().r() < firstr ||
3513 barrelcylinder == startlayer) {
3514 upstreamlayers.emplace_back(barrelcylinder, (
Layer*)
nullptr);
3517 if (barrelcylinder != startlayer &&
3518 (barrelcylinder->surfaceRepresentation().bounds().r() > firstr2 ||
3519 barrelcylinder == startlayer2)) {
3520 layers.emplace_back(barrelcylinder, (
Layer*)
nullptr);
3529 std::sort(upstreamlayers.begin(), upstreamlayers.end(),
GXF::LayerSort());
3533 const EventContext& ctx,
3581 addMaterial(ctx, cache, trajectory, refpar2, matEffects);
3597 bool hasmat =
false;
3598 int indexoffset = 0, lastmatindex = 0;
3599 std::vector<std::unique_ptr<GXFTrackState>> & oldstates = trajectory.
trackStates();
3613 for (
int i = 0;
i < (
int) oldstates.size();
i++) {
3614 if (oldstates[
i]->materialEffects() !=
nullptr) {
3623 if (firstsistate ==
nullptr) {
3624 if (oldstates[
i]->trackParameters() ==
nullptr) {
3625 std::unique_ptr<const TrackParameters> tmppar(
m_propagator->propagateParameters(
3628 oldstates[
i]->associatedSurface(),
3635 if (tmppar ==
nullptr)
return;
3637 oldstates[
i]->setTrackParameters(std::move(tmppar));
3639 firstsistate = oldstates[
i].get();
3641 lastsistate = oldstates[
i].get();
3649 if (lastsistate ==
nullptr) {
3650 throw std::logic_error(
"No track state");
3660 std::unique_ptr<const TrackParameters> tmppar(
m_propagator->propagateParameters(
3669 if (tmppar ==
nullptr)
return;
3681 indexoffset = lastmatindex;
3696 std::vector<std::pair<const Layer *, const Layer *>>
layers;
3697 std::vector<std::pair<const Layer *, const Layer *>> & upstreamlayers = trajectory.
upstreamMaterialLayers();
3711 const EventContext& ctx,
3717 if (refpar2 ==
nullptr) {
3727 std::vector<std::unique_ptr<GXFTrackState>> &
states = trajectory.
trackStates();
3728 std::vector<std::unique_ptr<GXFTrackState>> matstates;
3729 std::unique_ptr< const std::vector < const TrackStateOnSurface *>,
3730 void (*)(
const std::vector<const TrackStateOnSurface *> *) >
3731 matvec(
nullptr,&Trk::GlobalChi2Fitter::Cache::objVectorDeleter<TrackStateOnSurface>);
3732 bool matvec_used=
false;
3733 std::unique_ptr<TrackParameters> startmatpar1;
3734 std::unique_ptr<TrackParameters> startmatpar2;
3745 int npseudomuon1 = 0;
3746 int npseudomuon2 = 0;
3748 for (
auto & state :
states) {
3754 if (firstidhit ==
nullptr) {
3763 if (firsthit ==
nullptr) {
3764 firsthit = state->measurement();
3766 if (
tp ==
nullptr) {
3776 if (
tp ==
nullptr) {
3780 state->setTrackParameters(std::unique_ptr<const TrackParameters>(
tp));
3786 lasthit = state->measurement();
3792 if (firstidhit ==
nullptr) {
3793 firstidhit = state->measurement();
3796 if ((firstidpar ==
nullptr) && (
tp !=
nullptr)) {
3800 lastidhit = state->measurement();
3801 if (
tp !=
nullptr) {
3806 if (firstsiliconpar ==
nullptr) {
3807 firstsiliconpar =
tp;
3809 lastsiliconpar =
tp;
3821 if (firstmuonhit ==
nullptr) {
3822 firstmuonhit = state->measurement();
3823 if (
tp !=
nullptr) {
3827 lastmuonhit = state->measurement();
3828 if (
tp !=
nullptr) {
3834 if (meff->
deltaE() == 0) {
3835 if (firstcalopar ==
nullptr) {
3836 firstcalopar = state->trackParameters();
3838 lastcalopar = state->trackParameters();
3840 if (firstmatpar ==
nullptr) {
3841 firstmatpar = state->trackParameters();
3846 std::unique_ptr<TrackParameters> refpar;
3847 AmgVector(5) newpars = refpar2->parameters();
3854 newpars[0], newpars[1], newpars[2], newpars[3], newpars[4], std::nullopt
3857 if (firstmatpar !=
nullptr) {
3862 if ((startmatpar1 ==
nullptr) || ((firstidhit !=
nullptr) && (firstmuonhit !=
nullptr))) {
3868 const AmgVector(5) & newpars = startmatpar2->parameters();
3872 startmatpar2 = startmatpar2->associatedSurface().createUniqueTrackParameters(
3873 newpars[0], newpars[1], newpars[2], newpars[3],
3874 sign / std::sqrt(oldp * oldp + 2 * 100 *
MeV * std::sqrt(oldp * oldp +
mass *
mass) + 100 *
MeV * 100 *
MeV),
3879 AmgVector(5) newpars = startmatpar1->parameters();
3882 startmatpar1 = startmatpar1->associatedSurface().createUniqueTrackParameters(
3883 newpars[0], newpars[1], newpars[2], newpars[3], newpars[4], std::nullopt
3886 newpars = startmatpar2->parameters();
3889 startmatpar2 = startmatpar2->associatedSurface().createUniqueTrackParameters(
3890 newpars[0], newpars[1], newpars[2], newpars[3], newpars[4], std::nullopt
3898 refpar->momentum().unit()
3904 ATH_MSG_DEBUG(
"Obtaining upstream layers from Extrapolator");
3907 std::unique_ptr<const TrackParameters> tmppar;
3909 if (firstmuonhit !=
nullptr) {
3929 if (tmppar !=
nullptr) {
3930 destsurf = &tmppar->associatedSurface();
3935 if (matvec_used) cache.
m_matTempStore.push_back( std::move(matvec) );
3940 false, matEffects) );
3943 if (matvec && !matvec->empty()) {
3944 for (
int i = (
int)matvec->size() - 1;
i > -1;
i--) {
3949 std::unique_ptr<GXFMaterialEffects> meff = std::make_unique<GXFMaterialEffects>(*meot);
3950 const TrackParameters * newpars = (*matvec)[
i]->trackParameters() !=
nullptr ? (*matvec)[
i]->trackParameters()->
clone() :
nullptr;
3951 meff->setSigmaDeltaE(0);
3952 matstates.push_back(std::make_unique<GXFTrackState>(
3954 std::unique_ptr<const TrackParameters>(newpars)
3967 refpar->momentum().unit()
3973 ATH_MSG_DEBUG(
"Obtaining downstream ID layers from Extrapolator");
3975 std::unique_ptr<const TrackParameters> tmppar;
3976 std::unique_ptr<Surface> calosurf;
3977 if (firstmuonhit !=
nullptr) {
3998 if (tmppar !=
nullptr) {
4002 cylcalosurf =
static_cast<const CylinderSurface *
>(&tmppar->associatedSurface());
4007 disccalosurf =
static_cast<const DiscSurface *
>(&tmppar->associatedSurface());
4009 if (cylcalosurf !=
nullptr) {
4012 double radius = cylbounds.
r();
4014 calosurf = std::make_unique<CylinderSurface>(trans,
radius - 1, hlength);
4015 }
else if (disccalosurf !=
nullptr) {
4017 disccalosurf->
center().z() > 0 ?
4018 disccalosurf->
center().z() - 1 :
4019 disccalosurf->
center().z() + 1
4023 disccalosurf->
center().x(),
4024 disccalosurf->
center().y(),
4029 trans.translation() << newpos;
4032 double rmin = discbounds->
rMin();
4033 double rmax = discbounds->
rMax();
4034 calosurf = std::make_unique<DiscSurface>(trans, rmin, rmax);
4036 destsurf = calosurf.release();
4040 if (matvec_used) cache.
m_matTempStore.push_back( std::move(matvec) );
4042 ctx, *startmatpar2, *destsurf,
alongMomentum,
false, matEffects));
4043 matvec_used =
false;
4045 if (matvec && !matvec->empty()) {
4046 for (
const auto &
i : *matvec) {
4052 std::unique_ptr<GXFMaterialEffects> meff = std::make_unique<GXFMaterialEffects>(*meot);
4058 meff->setDeltaE(-5);
4061 meff->setScatteringSigmas(0, 0);
4064 meff->setSigmaDeltaE(50);
4067 const TrackParameters * newparams =
i->trackParameters() !=
nullptr ?
i->trackParameters()->clone() :
nullptr;
4069 matstates.push_back(std::make_unique<GXFTrackState>(
4071 std::unique_ptr<const TrackParameters>(newparams)
4083 if (cache.
m_calomat && (firstmuonhit !=
nullptr) && (firstidhit !=
nullptr)) {
4086 std::vector<MaterialEffectsOnTrack> calomeots =
m_calotool->extrapolationSurfacesAndEffects(
4095 if (calomeots.empty()) {
4100 std::unique_ptr<const TrackParameters> prevtrackpars =
unique_clone(lastidpar);
4101 if (lasthit == lastmuonhit) {
4102 for (
int i = 0;
i < (
int) calomeots.size();
i++) {
4105 std::unique_ptr<const TrackParameters> layerpar(
m_propagator->propagateParameters(
4108 calomeots[
i].associatedSurface(),
4115 if (layerpar ==
nullptr) {
4120 std::unique_ptr<GXFMaterialEffects> meff = std::make_unique<GXFMaterialEffects>(calomeots[
i]);
4123 lastcalopar = layerpar.get();
4127 double qoverp = layerpar->parameters()[
Trk::qOverP];
4128 double qoverpbrem = 0;
4132 (firstmuonpar !=
nullptr) &&
4133 std::abs(firstmuonpar->parameters()[
Trk::qOverP]) > 1.e-9
4135 qoverpbrem = firstmuonpar->parameters()[
Trk::qOverP];
4137 double sign = (qoverp > 0) ? 1 : -1;
4138 qoverpbrem =
sign / (1 / std::abs(qoverp) - std::abs(calomeots[
i].energyLoss()->deltaE()));
4141 const AmgVector(5) & newpar = layerpar->parameters();
4143 layerpar = layerpar->associatedSurface().createUniqueTrackParameters(
4144 newpar[0], newpar[1], newpar[2], newpar[3], qoverpbrem, std::nullopt
4146 meff->setdelta_p(1000 * (qoverpbrem - qoverp));
4149 matstates.push_back(std::make_unique<GXFTrackState>(
4151 std::unique_ptr<const TrackParameters>(layerpar !=
nullptr ? layerpar->clone() :
nullptr)
4153 prevtrackpars = std::move(layerpar);
4158 firsthit == firstmuonhit &&
4162 for (
int i = 0;
i < (
int) calomeots.size();
i++) {
4164 std::unique_ptr<const TrackParameters> layerpar(
m_propagator->propagateParameters(
4167 calomeots[
i].associatedSurface(),
4174 if (layerpar ==
nullptr) {
4179 std::unique_ptr<GXFMaterialEffects> meff = std::make_unique<GXFMaterialEffects>(calomeots[
i]);
4188 double qoverpbrem = layerpar->parameters()[
Trk::qOverP];
4193 (lastmuonpar !=
nullptr) &&
4194 std::abs(lastmuonpar->parameters()[
Trk::qOverP]) > 1.e-9
4198 double sign = (qoverpbrem > 0) ? 1 : -1;
4199 qoverp =
sign / (1 / std::abs(qoverpbrem) + std::abs(calomeots[
i].energyLoss()->deltaE()));
4202 meff->setdelta_p(1000 * (qoverpbrem - qoverp));
4203 const AmgVector(5) & newpar = layerpar->parameters();
4205 prevtrackpars = layerpar->associatedSurface().createUniqueTrackParameters(
4206 newpar[0], newpar[1], newpar[2], newpar[3], qoverp, std::nullopt
4210 matstates.insert(matstates.begin(), std::make_unique<GXFTrackState>(std::move(meff), std::move(layerpar)));
4215 if (lasthit == lastmuonhit && cache.
m_extmat) {
4216 std::unique_ptr<const Trk::TrackParameters> muonpar1;
4218 if (lastcalopar !=
nullptr) {
4238 if (muonpar1 !=
nullptr) {
4246 rot.col(2) = trackdir;
4248 trans.linear().matrix() << rot;
4249 trans.translation() << muonpar1->
position() - .1 * trackdir;
4252 std::unique_ptr<const TrackParameters> curvlinpar(
m_extrapolator->extrapolateDirectly(
4260 if (curvlinpar !=
nullptr) {
4261 muonpar1 = std::move(curvlinpar);
4265 muonpar1 = std::unique_ptr<TrackParameters>(lastcalopar->
clone());
4268 muonpar1 = std::unique_ptr<TrackParameters>(refpar->clone());
4273 if (muonpar1 !=
nullptr) {
4283 0) and (firstmuonhit !=
nullptr)) {
4295 std::abs(distsol.
first()) < std::abs(distsol.
second()) ?
4302 if (firstmuonpar !=
nullptr) {
4303 AmgVector(5) newpars = firstmuonpar->parameters();
4310 newpars[0], newpars[1], newpars[2], newpars[3], newpars[4], std::nullopt
4313 std::unique_ptr<const TrackParameters> tmppar(
m_propagator->propagateParameters(
4323 if (tmppar !=
nullptr) {
4324 muonpar1 = std::move(tmppar);
4330 ATH_MSG_DEBUG(
"Obtaining downstream layers from Extrapolator");
4331 if (matvec_used) cache.
m_matTempStore.push_back( std::move(matvec) );
4334 states.back()->associatedSurface(),
4338 matvec_used =
false;
4345 if (matvec && !matvec->empty()) {
4346 for (
int j = 0; j < (
int) matvec->size(); j++) {
4352 std::unique_ptr<GXFMaterialEffects> meff = std::make_unique<GXFMaterialEffects>(*meot);
4357 std::abs(meff->deltaE()) > 25 &&
4358 std::abs((*matvec)[j]->trackParameters()->position().z()) < 13000
4363 const TrackParameters * newparams = (*matvec)[j]->trackParameters() !=
nullptr ? (*matvec)[j]->trackParameters()->
clone() :
nullptr;
4365 matstates.push_back(std::make_unique<GXFTrackState>(
4367 std::unique_ptr<const TrackParameters>(newparams)
4377 if (firsthit == firstmuonhit && cache.
m_extmat && (firstcalopar !=
nullptr)) {
4378 std::unique_ptr<const Trk::TrackParameters> muonpar1;
4399 if (muonpar1 !=
nullptr) {
4407 rot.col(2) = trackdir;
4409 trans.linear().matrix() << rot;
4410 trans.translation() << muonpar1->
position() - .1 * trackdir;
4413 std::unique_ptr<const TrackParameters> curvlinpar(
m_extrapolator->extrapolateDirectly(
4421 if (curvlinpar !=
nullptr) {
4422 muonpar1 = std::move(curvlinpar);
4426 muonpar1 = std::unique_ptr<const TrackParameters>(firstcalopar->
clone());
4432 if (muonpar1 !=
nullptr) {
4443 ATH_MSG_DEBUG(
"Collecting upstream muon material from extrapolator");
4444 if (matvec_used) cache.
m_matTempStore.push_back( std::move(matvec) );
4447 states[0]->associatedSurface(),
4451 matvec_used =
false;
4453 if (matvec && !matvec->empty()) {
4454 ATH_MSG_DEBUG(
"Retrieved " << matvec->size() <<
" material states");
4456 for (
int j = 0; j < (
int) matvec->size(); j++) {
4459 if (meb !=
nullptr) {
4464 std::unique_ptr<GXFMaterialEffects> meff = std::make_unique<GXFMaterialEffects>(*meot);
4469 std::abs(meff->deltaE()) > 25 &&
4470 std::abs((*matvec)[j]->trackParameters()->position().z()) < 13000
4476 (*matvec)[j]->trackParameters() !=
nullptr
4477 ? (*matvec)[j]->trackParameters()->
clone()
4480 matstates.insert(matstates.begin(), std::make_unique<GXFTrackState>(
4482 std::unique_ptr<const TrackParameters>(tmpparams)
4495 std::vector<std::unique_ptr<GXFTrackState>> & newstates =
states;
4496 std::vector<std::unique_ptr<GXFTrackState>> oldstates = std::move(newstates);
4499 newstates.reserve(oldstates.size() + matstates.size());
4502 int firstlayerno = -1;
4512 bool addlayer =
true;
4514 while (addlayer && layerno < (
int) matstates.size()) {
4516 const TrackParameters *layerpar = matstates[layerno]->trackParameters();
4518 DistanceSolution distsol = oldstates[
i]->associatedSurface().straightLineDistanceEstimate(
4531 double trackimpact = std::abs(-refpar->position().x() * sinphi + refpar->position().y() * cosphi);
4533 if (trackimpact > cylinderradius - 5 *
mm) {
4539 if (
i == (
int) oldstates.size() - 1) {
4548 ATH_MSG_DEBUG(
"Material state with excessive scattering, skipping it");
4559 if (firstlayerno < 0) {
4560 firstlayerno = layerno;
4565 if ((layerpar !=
nullptr) && matEffects !=
pion && matEffects !=
muon) {
4567 const Layer *lay =
nullptr;
4569 if (tvol !=
nullptr) {
4585 if (matvec_used) cache.
m_matTempStore.push_back( std::move(matvec) );
4610 std::unique_ptr<const TrackParameters> per(
m_extrapolator->extrapolate(
4611 Gaudi::Hive::currentContext(),param, tmppersf,
oppositeMomentum,
false, matEffects));
4613 if (per ==
nullptr) {
4614 ATH_MSG_DEBUG(
"Cannot make Perigee with starting parameters");
4618 if(std::abs(per->position().z())>5000.) {
4619 ATH_MSG_WARNING(
"Pathological perigee well outside of tracking detector!! Returning nullptr");
4627 const EventContext& ctx,
4634 ATH_MSG_DEBUG(
"--> entering GlobalChi2Fitter::myfit_helper");
4662 if (trajectory.
nDOF() < 0) {
4672 ATH_MSG_WARNING(
"Attempt to apply material corrections with q/p=0, reject track");
4686 std::unique_ptr<const TrackParameters> per =
makePerigee(cache, param, matEffects);
4702 ATH_MSG_WARNING(
"Tracking Geometry Service and/or Material Updator Tool not configured");
4713 addMaterial(ctx, cache, trajectory, per !=
nullptr ? per.get() : ¶m, matEffects);
4716 ctx, cache, trajectory, per !=
nullptr ? per.get() : ¶m, matEffects);
4728 for (std::vector<std::unique_ptr<GXFTrackState>>::
iterator it =
4735 (**it).materialEffects()->deltaE() == 0
4737 scatstate2 = (*it).get();
4742 scatstate = (*it).get();
4749 if ((scatstate ==
nullptr) || (scatstate2 ==
nullptr)) {
4750 throw std::logic_error(
"Invalid scatterer");
4757 trajectory.
trackStates()[nstates / 2 - 1]->position() +
4763 std::unique_ptr<const TrackParameters> nearestpar;
4764 double mindist = 99999;
4765 std::vector < GXFTrackState * >mymatvec;
4768 if ((*it).trackParameters() ==
nullptr) {
4774 (*it).trackParameters()->position(),
4775 (*it).trackParameters()->momentum().unit())
4784 (((*it).measurement() !=
nullptr) && insideid) || (
4785 ((*it).materialEffects() !=
nullptr) &&
4787 (*it).materialEffects()->deltaE() == 0 ||
4788 ((*it).materialEffects()->sigmaDeltaPhi() == 0 &&
4790 (*it).materialEffects()->deltaPhi() != 0
4794 double dist = ((*it).trackParameters()->position() -
vertex).
perp();
4795 if (dist < mindist) {
4803 if (((*it).materialEffects() !=
nullptr) &&
distance > 0) {
4804 mymatvec.push_back(
it.get());
4808 if (nearestpar ==
nullptr) {
4812 for (
auto & state : mymatvec) {
4814 const Surface &matsurf = state->associatedSurface();
4816 nearestpar->position(), nearestpar->momentum().unit());
4824 std::unique_ptr<const TrackParameters> tmppar(
m_propagator->propagateParameters(
4834 if (tmppar ==
nullptr) {
4846 if (tmppar ==
nullptr) {
4856 AmgVector(5) newpars = tmppar->parameters();
4858 if (state->materialEffects()->sigmaDeltaE() > 0) {
4859 newpars[
Trk::qOverP] += .001 * state->materialEffects()->delta_p();
4862 double de = std::abs(state->materialEffects()->deltaE());
4864 double newp2 = oldp * oldp - 2 * de * std::sqrt(
mass *
mass + oldp * oldp) + de * de;
4870 nearestpar = tmppar->associatedSurface().createUniqueTrackParameters(
4871 newpars[0], newpars[1], newpars[2], newpars[3], newpars[4], std::nullopt
4875 std::unique_ptr<Trk::TrackParameters> tmpPars(
m_propagator->propagateParameters(
4886 if (tmpPars !=
nullptr) {
4887 per.reset(
static_cast < const Perigee *
>(tmpPars.release()));
4892 double oldp = 1. / std::abs(per->parameters()[
Trk::qOverP]);
4894 double newp = std::sqrt(oldp * oldp + 2 * toteloss * std::sqrt(oldp * oldp +
mass *
mass) + toteloss * toteloss);
4898 per = per->associatedSurface().createUniqueTrackParameters(
4903 if (per ==
nullptr) {
4920 }
else if (per ==
nullptr) {
4938 per = per->associatedSurface().createUniqueTrackParameters(
4945 if (per !=
nullptr) {
4955 Eigen::MatrixXd a_inv;
4956 a.resize(nfitpar, nfitpar);
4961 derivPool.setZero();
4963 for (std::unique_ptr<GXFTrackState> & state : trajectory.
trackStates()) {
4964 if (state->materialEffects() !=
nullptr) {
4967 state->setDerivatives(derivPool);
4970 bool doderiv =
true;
5003 double redchi2 = (trajectory.
nDOF() > 0) ? trajectory.
chi2() / trajectory.
nDOF() : 0;
5004 double prevredchi2 = (trajectory.
nDOF() > 0) ? trajectory.
prevchi2() / trajectory.
nDOF() : 0;
5013 (redchi2 < prevredchi2 &&
5014 (redchi2 > prevredchi2 - 1 || redchi2 < 2)) ||
5015 nsihits + ntrthits == nhits
5020 if (
it != 1 || nsihits != 0 || trajectory.
nDOF() <= 0 || trajectory.
chi2() / trajectory.
nDOF() <= 3) {
5035 if (ntrtprechits+ntrttubehits) {
5036 phf =
float(ntrtprechits)/
float(ntrtprechits+ntrttubehits);
5038 if (phf<m_minphfcut && it>=3) {
5039 if ((ntrtprechits+ntrttubehits)>=15) {
5044 <<
" | nTRTPrecHits = " << ntrtprechits
5045 <<
" | nTRTTubeHits = " << ntrttubehits
5046 <<
" | nOutliers = "
5066 if (trajectory.
prefit() == 0) {
5070 Eigen::LLT < Eigen::MatrixXd > lltOfW(
a);
5071 if (lltOfW.info() == Eigen::Success) {
5075 int ncols =
a.cols();
5076 Amg::MatrixX weightInvAMG = Amg::MatrixX::Identity(ncols, ncols);
5077 a_inv = lltOfW.solve(weightInvAMG);
5097 if (finaltrajectory != &trajectory) {
5098 delete finaltrajectory;
5108 if (nperpars == 5) {
5109 for (
int i = 0;
i <
a.cols();
i++) {
5110 a_inv(4,
i) *= .001;
5111 a_inv(
i, 4) *= .001;
5115 int scatterPos = nperpars + 2 * nscat;
5116 for (
int bremno = 0; bremno < nbrem; bremno++, scatterPos++) {
5117 for (
int i = 0;
i <
a.cols();
i++) {
5118 a_inv(scatterPos,
i) *= .001;
5119 a_inv(
i, scatterPos) *= .001;
5126 for (
int i = 0;
i < nperparams;
i++) {
5127 for (
int j = 0; j < nperparams; j++) {
5128 (errmat) (j,
i) = a_inv(j,
i);
5133 (errmat) (4, 4) = 1
e-20;
5137 std::unique_ptr<const TrackParameters> measper(
5139 perpars[0], perpars[1], perpars[2], perpars[3], perpars[4], std::move(errmat)
5149 std::unique_ptr<Track>
track =
nullptr;
5151 if (finaltrajectory->
prefit() > 0) {
5152 if (finaltrajectory != &trajectory) {
5153 delete finaltrajectory;
5172 (
track !=
nullptr) && (
5173 track->fitQuality()->numberDoF() != 0 &&
5174 track->fitQuality()->chiSquared() /
track->fitQuality()->numberDoF() >
cut
5177 track.reset(
nullptr);
5181 if (
track ==
nullptr) {
5185 if (finaltrajectory != &trajectory) {
5186 delete finaltrajectory;
5189 return track.release();
5193 const EventContext& ctx,
5204 std::vector<std::unique_ptr<GXFTrackState>> &
states = trajectory.
trackStates();
5219 int nmeas = (
int)
res.size();
5222 bool scatwasupdated =
false;
5225 int bremno_maxbrempull = 0;
5226 double maxbrempull = 0;
5228 for (
int hitno = 0; hitno < (
int)
states.size(); hitno++) {
5229 std::unique_ptr<GXFTrackState> & state =
states[hitno];
5239 trajectory.
nDOF() != 0 &&
5240 std::abs((trajectory.
prevchi2() - trajectory.
chi2()) / trajectory.
nDOF()) < 15 &&
5241 !state->associatedSurface().isFree() &&
5243 (nperpars == 0 || nidhits > 0) &&
5244 (!state->isRecalibrated())
5249 std::unique_ptr<const PseudoMeasurementOnTrack> newpseudo = std::make_unique<const PseudoMeasurementOnTrack>(
5255 state->setMeasurement(std::move(newpseudo));
5256 measbase = state->measurement();
5259 double *
errors = state->measurementErrors();
5263 for (
int i = 0;
i < 5;
i++) {
5276 res[measno] = residuals[
i];
5278 if (
i == 2 && std::abs(std::abs(
res[measno]) - 2 *
M_PI) < std::abs(
res[measno])) {
5279 if (
res[measno] < 0) {
5288 double *
errors = state->measurementErrors();
5289 for (
int i = 0;
i < 5;
i++) {
5299 ((trajectory.
prefit() == 0) || state->materialEffects()->deltaE() == 0)
5301 double deltaphi = state->materialEffects()->deltaPhi();
5302 double measdeltaphi = state->materialEffects()->measuredDeltaPhi();
5303 double sigmadeltaphi = state->materialEffects()->sigmaDeltaPhi();
5304 double deltatheta = state->materialEffects()->deltaTheta();
5305 double sigmadeltatheta = state->materialEffects()->sigmaDeltaTheta();
5307 if (trajectory.
prefit() != 1) {
5308 b[nperpars + 2 * scatno] -= (deltaphi - measdeltaphi) / (sigmadeltaphi * sigmadeltaphi);
5309 b[nperpars + 2 * scatno + 1] -= deltatheta / (sigmadeltatheta * sigmadeltatheta);
5311 b[nperpars + scatno] -= deltatheta / (sigmadeltatheta * sigmadeltatheta);
5315 deltaphi * deltaphi / (sigmadeltaphi * sigmadeltaphi) +
5316 deltatheta * deltatheta / (sigmadeltatheta * sigmadeltatheta)
5322 if ((state->materialEffects() !=
nullptr) && state->materialEffects()->sigmaDeltaE() > 0) {
5323 double averagenergyloss = std::abs(state->materialEffects()->deltaE());
5325 const double qoverp = limitInversePValue(qoverpbrem - state->materialEffects()->delta_p());
5326 const double pbrem = 1. / std::abs(qoverpbrem);
5327 const double p = 1. / std::abs(qoverp);
5328 const double mass = .001 * trajectory.
mass();
5330 const double bremEnergy = std::sqrt(pbrem * pbrem +
mass *
mass);
5332 res[nmeas - nbrem + bremno] = .001 * averagenergyloss -
energy + bremEnergy;
5334 double sigde = state->materialEffects()->sigmaDeltaE();
5335 double sigdepos = state->materialEffects()->sigmaDeltaEPos();
5336 double sigdeneg = state->materialEffects()->sigmaDeltaENeg();
5338 error[nmeas - nbrem + bremno] = .001 * state->materialEffects()->sigmaDeltaE();
5340 if (state->materialEffects()->isKink()) {
5341 maxbrempull = -999999999;
5342 state_maxbrempull =
nullptr;
5348 (trajectory.
prefit() == 0) &&
5350 sigde != sigdepos &&
5353 double elosspull =
res[nmeas - nbrem + bremno] / (.001 * sigde);
5355 if (trajectory.
mass() > 100) {
5356 if (elosspull < -1) {
5357 state->materialEffects()->setSigmaDeltaE(sigdepos);
5358 }
else if (elosspull > 1) {
5359 state->materialEffects()->setSigmaDeltaE(sigdeneg);
5362 error[nmeas - nbrem + bremno] = .001 * state->materialEffects()->sigmaDeltaE();
5365 !state->materialEffects()->isKink() && (
5366 (elosspull < -.2 &&
m_fixbrem == -1 && elosspull < maxbrempull) ||
5370 bremno_maxbrempull = bremno;
5371 state_maxbrempull = state.get();
5372 maxbrempull = elosspull;
5381 (trajectory.
prefit() == 0) &&
5382 state->materialEffects()->sigmaDeltaPhi() == 0 &&
5383 state->materialEffects()->isMeasuredEloss() &&
5384 res[nmeas - nbrem + bremno] / (.001 * state->materialEffects()->sigmaDeltaEAve()) > 2.5
5388 :
states[hitno - 2]->trackParameters();
5391 std::vector<MaterialEffectsOnTrack> calomeots =
5400 if (calomeots.size() == 3) {
5401 averagenergyloss = std::abs(calomeots[1].energyLoss()->deltaE());
5402 double newres = .001 * averagenergyloss -
energy + bremEnergy;
5403 double newerr = .001 * calomeots[1].energyLoss()->sigmaDeltaE();
5405 if (std::abs(newres / newerr) < std::abs(
res[nmeas - nbrem + bremno] /
error[nmeas - nbrem + bremno])) {
5406 ATH_MSG_DEBUG(
"Changing from measured to parametrized energy loss");
5408 state->materialEffects()->setEloss(std::unique_ptr<EnergyLoss>(calomeots[1].energyLoss()->
clone()));
5409 state->materialEffects()->setSigmaDeltaE(calomeots[1].energyLoss()->sigmaDeltaE());
5410 res[nmeas - nbrem + bremno] = newres;
5411 error[nmeas - nbrem + bremno] = newerr;
5415 state->materialEffects()->setMeasuredEloss(
false);
5423 for (; measno < nmeas; measno++) {
5424 if (
error[measno] == 0) {
5431 if (!doderiv && (scatwasupdated)) {
5435 double oldchi2 = trajectory.
chi2();
5439 double oldredchi2 = (trajectory.
nDOF() > 0) ? oldchi2 / trajectory.
nDOF() : 0;
5440 double newredchi2 = (trajectory.
nDOF() > 0) ?
chi2 / trajectory.
nDOF() : 0;
5443 trajectory.
prefit() > 0 && (
5444 (newredchi2 < 2 &&
it != 0) ||
5445 (newredchi2 < oldredchi2 + .1 && std::abs(newredchi2 - oldredchi2) < 1 &&
it != 1)
5451 double maxdiff = (nsihits != 0 && nsihits + ntrthits == nhits &&
chi2 < oldchi2) ? 200 : 1.;
5453 int miniter = (nsihits != 0 && nsihits + ntrthits == nhits) ? 1 : 2;
5463 if ((state_maxbrempull !=
nullptr) && trajectory.
converged()) {
5471 double olderror =
error[nmeas - nbrem + bremno_maxbrempull];
5475 if (
a.cols() != nfitpars) {
5479 for (
int i = 0;
i < nfitpars;
i++) {
5480 if (weightderiv(nmeas - nbrem + bremno_maxbrempull,
i) == 0) {
5484 for (
int j =
i; j < nfitpars; j++) {
5488 weightderiv(nmeas - nbrem + bremno_maxbrempull,
i) *
5489 weightderiv(nmeas - nbrem + bremno_maxbrempull, j) *
5490 (1 - olderror * olderror / (newerror * newerror))
5494 weightderiv(nmeas - nbrem + bremno_maxbrempull,
i) *= olderror / newerror;
5508 std::vector<std::unique_ptr<GXFTrackState>> &
states = trajectory.
trackStates();
5521 int nmeas = (
int) weightderiv.rows();
5525 for (std::unique_ptr<GXFTrackState> & state :
states) {
5528 ((state->materialEffects() ==
nullptr) || state->materialEffects()->sigmaDeltaE() <= 0)
5535 const auto [scatmin, scatmax] = std::minmax(scatno, nscatupstream);
5536 const auto [bremmin, bremmax] = std::minmax(bremno, nbremupstream);
5539 double sinstereo = 0;
5542 sinstereo = state->sinStereo();
5545 double cosstereo = (sinstereo == 0) ? 1. : std::sqrt(1 - sinstereo * sinstereo);
5547 for (
int i = 0;
i < 5;
i++) {
5559 weightderiv.row(measno).head(
cols) =
5560 (derivatives.row(0).head(
cols) * cosstereo +
5561 sinstereo * derivatives.row(1).head(
cols)) /
5564 weightderiv.row(measno).head(
cols) = derivatives.row(
i).head(
cols) /
error[measno];
5568 for (
int j = scatmin; j < scatmax; j++) {
5569 int index = nperparams + ((trajectory.
prefit() != 1) ? 2 * j : j);
5570 double thisderiv = 0;
5574 if (
i == 0 && sinstereo != 0) {
5575 thisderiv =
sign * (derivatives(0,
index) * cosstereo + sinstereo * derivatives(1,
index));
5580 weightderiv(measno,
index) = thisderiv /
error[measno];
5582 if (trajectory.
prefit() != 1) {
5585 if (
i == 0 && sinstereo != 0) {
5586 thisderiv =
sign * (derivatives(0,
index) * cosstereo + sinstereo * derivatives(1,
index));
5591 weightderiv(measno,
index) = thisderiv /
error[measno];
5595 for (
int j = bremmin; j < bremmax; j++) {
5596 double thisderiv = 0;
5597 int index = j + nperparams + 2 * nscat;
5598 if (
i == 0 && sinstereo != 0) {
5599 thisderiv = derivatives(0,
index) * cosstereo + sinstereo * derivatives(1,
index);
5601 thisderiv = derivatives(
i,
index);
5603 weightderiv(measno,
index) = thisderiv /
error[measno];
5609 double *
errors = state->measurementErrors();
5610 for (
int i = 0;
i < 5;
i++) {
5617 ((trajectory.
prefit() == 0) || state->materialEffects()->deltaE() == 0)
5622 if ((state->materialEffects() !=
nullptr) && state->materialEffects()->sigmaDeltaE() > 0) {
5624 double qoverpbrem = limitInversePValue(1000 * state->trackParameters()->parameters()[
Trk::qOverP]);
5625 double qoverp = limitInversePValue(qoverpbrem - state->materialEffects()->delta_p());
5627 double mass = .001 * trajectory.
mass();
5629 const auto thisMeasurementIdx{nmeas - nbrem + bremno};
5633 auto multiplier = [] (
double mass,
double qOverP){
5636 const auto qoverpTerm {multiplier(
mass, qoverp) /
error[thisMeasurementIdx]};
5637 const auto qoverpBremTerm {multiplier(
mass, qoverpbrem) /
error[thisMeasurementIdx]};
5640 weightderiv(thisMeasurementIdx, 4) = qoverpBremTerm - qoverpTerm;
5643 const auto bremNoBase = nperparams + 2 * nscat;
5644 if (bremno < nbremupstream) {
5645 weightderiv(thisMeasurementIdx, bremNoBase + bremno) = qoverpTerm;
5646 for (
int bremno2 = bremno + 1; bremno2 < nbremupstream; bremno2++) {
5647 weightderiv(thisMeasurementIdx, bremNoBase + bremno2) = qoverpTerm - qoverpBremTerm;
5650 weightderiv(thisMeasurementIdx, bremNoBase + bremno) = qoverpBremTerm;
5651 for (
int bremno2 = nbremupstream; bremno2 < bremno; bremno2++) {
5652 weightderiv(thisMeasurementIdx, bremNoBase + bremno2) = qoverpBremTerm - qoverpTerm;
5661 const EventContext& ctx,
5676 double oldchi2 = trajectory.
chi2();
5677 double oldredchi2 = (trajectory.
nDOF() > 0) ? oldchi2 / trajectory.
nDOF() : 0;
5696 double newredchi2 = (trajectory.
nDOF() > 0) ? trajectory.
chi2() / trajectory.
nDOF() : 0;
5699 "=" << oldredchi2 <<
" new chi2: " << trajectory.
chi2() <<
"/" <<
5700 trajectory.
nDOF() <<
"=" << newredchi2);
5708 std::vector < std::pair < double, double >>&scatsigmas = trajectory.
scatteringSigmas();
5710 int nmeas = (
int)
res.size();
5722 for (
int i = 0;
i < nperpars;
i++) {
5730 std::unique_ptr<GXFTrackState> & state = trajectory.
trackStates()[
i];
5732 if (meff ==
nullptr) {
5733 measno += state->numberOfMeasuredParameters();
5735 if (meff !=
nullptr) {
5737 && ((trajectory.
prefit() == 0) || meff->
deltaE() == 0)) {
5738 int scatterPos = nperpars + 2 * scatno;
5739 if (
i < nupstreamstates) {
5753 if (
i < nupstreamstates) {
5768 if (
a.cols() != nfitpars) {
5772 for (
int k = 0;
k < nfitpars;
k++) {
5774 int maxmeas = nmeas - nbrem;
5778 for (measno = minmeas; measno < maxmeas; measno++) {
5780 res[measno] * (1. /
error[measno]) * weight_deriv(measno,
k);
5784 if (
k == 4 ||
k >= nperpars + scatpars) {
5785 for (measno = nmeas - nbrem; measno < nmeas; measno++) {
5786 b[
k] +=
res[measno] * (1. /
error[measno]) * weight_deriv(measno,
k);
5791 for (
int l =
k;
l < nfitpars;
l++) {
5798 for (measno = minmeas; measno < maxmeas; measno++) {
5799 tmp += weight_deriv(measno,
k) * weight_deriv(measno,
l);
5801 a.fillSymmetric(
l,
k,
tmp);
5809 for (
int k = nperpars;
k < nperpars + scatpars;
k += 2) {
5810 a(
k,
k) += 1. / (scatsigmas[scatno].first * scatsigmas[scatno].first);
5811 a(
k + 1,
k + 1) += 1. / (scatsigmas[scatno].second * scatsigmas[scatno].second);
5815 for (
int measno = nmeas - nbrem; measno < nmeas; measno++) {
5816 for (
int k = 4;
k < nfitpars;
k++) {
5818 k = nperpars + scatpars;
5821 for (
int l =
k;
l < nfitpars;
l++) {
5823 l = nperpars + scatpars;
5825 double tmp =
a(
l,
k) + weight_deriv(measno,
k) * weight_deriv(measno,
l);
5826 a.fillSymmetric(
l,
k,
tmp);
5832 unsigned int scatno = 0;
5833 bool weightchanged =
false;
5835 for (std::unique_ptr<GXFTrackState> & thisstate : trajectory.
trackStates()) {
5838 if (meff !=
nullptr) {
5842 plsurf =
static_cast < const PlaneSurface *
>(&thisstate->associatedSurface());
5843 if (meff->
deltaE() == 0 || ((trajectory.
prefit() == 0) && (plsurf !=
nullptr))) {
5844 weightchanged =
true;
5846 if (
a.cols() != nfitpars) {
5850 int scatNoIndex = 2 * scatno + nperpars;
5852 if (trajectory.
prefit() == 0) {
5853 if (thisstate->materialEffects()->sigmaDeltaPhi() != 0) {
5856 message <<
"scatno is out of range " << scatno <<
" !< " << cache.
m_phiweight.size();
5857 throw std::range_error(
message.str());
5861 a(scatNoIndex, scatNoIndex) /= cache.
m_phiweight[scatno];
5866 }
else if (
it == 1) {
5868 }
else if (
it <= 3) {
5870 }
else if (
it <= 6) {
5876 a(scatNoIndex, scatNoIndex) *= cache.
m_phiweight[scatno];
5880 else if (trajectory.
prefit() >= 2) {
5881 if (newredchi2 > oldredchi2 - 1 && newredchi2 < oldredchi2) {
5882 a(scatNoIndex, scatNoIndex) *= 1.0001;
5883 a(scatNoIndex + 1, scatNoIndex + 1) *= 1.0001;
5884 }
else if (newredchi2 > oldredchi2 - 25 && newredchi2 < oldredchi2) {
5885 a(scatNoIndex, scatNoIndex) *= 1.001;
5886 a(scatNoIndex + 1, scatNoIndex + 1) *= 1.0001;
5888 a(scatNoIndex, scatNoIndex) *= 1.1;
5889 a(scatNoIndex + 1, scatNoIndex + 1) *= 1.001;
5895 thisstate->materialEffects()->sigmaDeltaPhi() != 0 &&
5896 ((trajectory.
prefit() == 0) || thisstate->materialEffects()->deltaE() == 0)
5904 (trajectory.
prefit() == 2) &&
5907 (newredchi2 < oldredchi2 - 25 || newredchi2 > oldredchi2)
5912 if (doderiv || weightchanged) {
5917 if ((trajectory.
prefit() == 0) && nsihits + ntrthits != nhits) {
5918 unsigned int scatno = 0;
5920 if (
a.cols() != nfitpars) {
5924 for (std::unique_ptr<GXFTrackState> & thisstate : trajectory.
trackStates()) {
5925 if ((thisstate->materialEffects() !=
nullptr) && thisstate->materialEffects()->sigmaDeltaPhi() != 0) {
5928 message <<
"scatno is out of range " << scatno <<
" !< " << cache.
m_phiweight.size();
5929 throw std::range_error(
message.str());
5935 plsurf =
static_cast<const PlaneSurface *
>(&thisstate->associatedSurface());
5937 if (thisstate->materialEffects()->deltaE() == 0 || (plsurf !=
nullptr)) {
5938 int scatNoIndex = 2 * scatno + nperpars;
5939 a(scatNoIndex, scatNoIndex) /= cache.
m_phiweight[scatno];
5943 if (thisstate->materialEffects()->sigmaDeltaPhi() != 0) {
5956 (newredchi2 < 2 || (newredchi2 < oldredchi2 && newredchi2 > oldredchi2 - .5)) &&
5957 (trajectory.
prefit() == 0)
5973 double d0 = refpar->parameters()[
Trk::d0];
5974 double z0 = refpar->parameters()[
Trk::z0];
5977 double qoverp = refpar->parameters()[
Trk::qOverP];
5982 Eigen::LLT<Eigen::MatrixXd> llt(lu_m);
5985 if (llt.info() == Eigen::Success) {
6005 std::vector < std::pair < double, double >>&scatangles = trajectory.
scatteringAngles();
6006 std::vector < double >&delta_ps = trajectory.
brems();
6008 for (
int i = 0;
i < nscat;
i++) {
6009 scatangles[
i].first +=
result[2 *
i + nperparams];
6010 scatangles[
i].second +=
result[2 *
i + nperparams + 1];
6013 for (
int i = 0;
i < nbrem;
i++) {
6014 delta_ps[
i] +=
result[nperparams + 2 * nscat +
i];
6017 std::unique_ptr<const TrackParameters> newper(
6034 const EventContext& evtctx
6045 if (!splitProbContainer.
isValid()) {
6049 std::vector<std::unique_ptr<GXFTrackState>> &
states = trajectory.
trackStates();
6056 for (
size_t stateno = 0; stateno <
states.size(); stateno++) {
6061 measno +=
states[stateno-1]->numberOfMeasuredParameters();
6064 std::unique_ptr<GXFTrackState> & state =
states[stateno];
6076 const auto *
const rot =
static_cast<const RIO_OnTrack *
>(pMeas);
6088 if (!splitProb.isSplit()) {
6089 ATH_MSG_DEBUG(
"Pixel cluster is not split so no need to update" );
6093 std::unique_ptr < const RIO_OnTrack > newrot;
6094 double *olderror = state->measurementErrors();
6097 double newerror[5] = {-1,-1,-1,-1,-1};
6098 double newres[2] = {-1,-1};
6100 newrot.reset(
m_ROTcreator->correct(*prd, *trackpars, evtctx));
6107 newerror[0] = std::sqrt(
covmat(0, 0));
6108 newres[0] = newrot->localParameters()[
Trk::locX] - trackpars->parameters()[
Trk::locX];
6109 newerror[1] = std::sqrt(
covmat(1, 1));
6110 newres[1] = newrot->localParameters()[
Trk::locY] - trackpars->parameters()[
Trk::locY];
6112 if (
a.cols() != nfitpars) {
6117 for(
int k =0;
k<2;
k++ ){
6118 double oldres =
res[measno+
k];
6119 res[measno+
k] = newres[
k];
6120 err[measno+
k] = newerror[
k];
6122 for (
int i = 0;
i < nfitpars;
i++) {
6123 if (weightderiv(measno+
k,
i) == 0) {
6127 b[
i] -= weightderiv(measno+
k,
i) * (oldres / olderror[
k] - (newres[
k] * olderror[
k]) / (newerror[
k] * newerror[
k]));
6129 for (
int j =
i; j < nfitpars; j++) {
6133 weightderiv(measno+
k,
i) *
6134 weightderiv(measno+
k, j) *
6135 ((olderror[
k] * olderror[
k]) / (newerror[
k] * newerror[
k]) - 1)
6139 weightderiv(measno+
k,
i) *= olderror[
k] / newerror[
k];
6143 state->setMeasurement(std::move(newrot));
6144 state->setMeasurementErrors(newerror);
6159 const EventContext& ctx
6167 std::vector<std::unique_ptr<GXFTrackState>> &
states = trajectory.
trackStates();
6173 if (
a.cols() != nfitpars) {
6181 bool outlierremoved =
false;
6182 bool hitrecalibrated =
false;
6184 for (
int stateno = 0; stateno < (
int)
states.size(); stateno++) {
6185 std::unique_ptr<GXFTrackState> & state =
states[stateno];
6193 std::abs(state->trackParameters()->parameters()[
Trk::driftRadius]) > 1.05 * state->associatedSurface().bounds().r()
6198 outlierremoved =
true;
6200 double *
errors = state->measurementErrors();
6201 double olderror =
errors[0];
6205 for (
int i = 0;
i < nfitpars;
i++) {
6206 if (weightderiv(measno,
i) == 0) {
6210 b[
i] -=
res[measno] * weightderiv(measno,
i) / olderror;
6212 for (
int j =
i; j < nfitpars; j++) {
6215 a(
i, j) - weightderiv(measno,
i) * weightderiv(measno, j)
6218 weightderiv(measno,
i) = 0;
6222 }
else if (trtrecal) {
6223 double *
errors = state->measurementErrors();
6224 double olderror =
errors[0];
6226 const auto *
const thisMeasurement{state->measurement()};
6232 if (oldrot->prepRawData() !=
nullptr) {
6233 double dcradius = oldrot->prepRawData()->localPosition()[
Trk::driftRadius];
6235 double trackradius = state->trackParameters()->parameters()[
Trk::driftRadius];
6237 std::unique_ptr<const Trk::RIO_OnTrack> newrot =
nullptr;
6238 double distance = std::abs(std::abs(trackradius) - dcradius);
6240 if (
distance < scalefactor * dcerror && (olderror > 1. || trackradius * oldradius < 0)) {
6241 newrot.reset(
m_ROTcreator->correct(*oldrot->prepRawData(), *state->trackParameters(), ctx));
6242 }
else if (
distance > scalefactor * dcerror && olderror < 1.) {
6243 newrot.reset(
m_broadROTcreator->correct(*oldrot->prepRawData(), *state->trackParameters(), ctx));
6246 if (newrot !=
nullptr) {
6248 hitrecalibrated =
true;
6252 if ((measno < 0) or (measno >= (
int)
res.size())) {
6253 throw std::runtime_error(
6254 "'res' array index out of range in TrkGlobalChi2Fitter/src/GlobalChi2Fitter.cxx:" +
std::to_string(__LINE__)
6258 double oldres =
res[measno];
6259 double newres = newradius - state->trackParameters()->parameters()[
Trk::driftRadius];
6261 state->setMeasurement(std::move(newrot));
6265 for (
int i = 0;
i < nfitpars;
i++) {
6266 if (weightderiv(measno,
i) == 0) {
6270 b[
i] -= weightderiv(measno,
i) * (oldres / olderror - (newres * olderror) / (newerror * newerror));
6272 for (
int j =
i; j < nfitpars; j++) {
6279 i < nperpars + 2 * nscats &&
6280 (
i - nperpars) % 2 == 0
6287 a(
i, j) + weightderiv(measno,
i) * weightderiv(measno, j) * ((olderror * olderror) / (newerror * newerror) - 1) *
weight
6290 weightderiv(measno,
i) *= olderror / newerror;
6293 res[measno] = newres;
6294 err[measno] = newerror;
6303 measno += state->numberOfMeasuredParameters();
6307 if (trajectory.
nDOF() < 0) {
6312 if (outlierremoved || hitrecalibrated) {
6321 const EventContext& ctx,
6329 bool trackok =
false;
6331 std::unique_ptr < GXFTrajectory > cleanup_oldtrajectory;
6333 std::unique_ptr < GXFTrajectory > cleanup_newtrajectory;
6339 while (!trackok && oldtrajectory->
nDOF() > 0) {
6341 std::vector<std::unique_ptr<GXFTrackState>> &
states = oldtrajectory->
trackStates();
6349 if (nhits != nsihits) {
6353 double maxsipull = -1;
6355 int hitno_maxsipull = -1;
6356 int measno_maxsipull = -1;
6357 int stateno_maxsipull = 0;
6369 for (
int stateno = 0; stateno < (
int)
states.size(); stateno++) {
6370 std::unique_ptr<GXFTrackState> & state =
states[stateno];
6376 double *
errors = state->measurementErrors();
6378 const Amg::MatrixX & hitcov = state->measurement()->localCovariance();
6379 double sinstereo = state->sinStereo();
6380 double cosstereo = (sinstereo == 0) ? 1 : std::sqrt(1 - sinstereo * sinstereo);
6381 double weight1 = -1;
6383 if (hitcov(0, 0) > trackcov(0, 0)) {
6384 if (sinstereo == 0) {
6388 trackcov(0, 0) * cosstereo * cosstereo + 2 *
6389 trackcov(1, 0) * cosstereo * sinstereo + trackcov(1, 1) * sinstereo * sinstereo
6400 double sipull1 = weight1 > 0 ? std::abs(
res[measno] / std::sqrt(weight1)) : -1;
6403 std::abs(
res[measno + 1] / std::sqrt(weight2)) :
6406 sipull1 =
std::max(sipull1, sipull2);
6408 if (sipull1 > maxsipull) {
6409 maxsipull = sipull1;
6410 measno_maxsipull = measno;
6411 state_maxsipull = state.get();
6412 stateno_maxsipull = stateno;
6413 hitno_maxsipull = hitno;
6424 measno += state->numberOfMeasuredParameters();
6428 double maxpull = maxsipull;
6430 ATH_MSG_DEBUG(
" maxsipull: " << maxsipull <<
" hitno_maxsipull: " <<
6431 hitno_maxsipull <<
" n3sigma: " << n3sigma <<
" cut: " <<
cut <<
" cut2: " << cut2);
6442 state_maxsipull = oldtrajectory->
trackStates()[stateno_maxsipull].get();
6445 const auto *
const rot =
static_cast<const RIO_OnTrack *
>(pMeas);
6448 std::unique_ptr < const RIO_OnTrack > broadrot;
6453 Amg::VectorX parameterVector = trackpar_maxsipull->parameters();
6454 std::unique_ptr<const TrackParameters> trackparForCorrect(
6463 state_maxsipull->trackCovariance())
6467 newerror[0] = newerror[1] = newerror[2] = newerror[3] = newerror[4] = -1;
6468 double newpull = -1;
6469 double newpull1 = -1;
6470 double newpull2 = -1;
6471 double newres1 = -1;
6472 double newres2 = -1;
6473 double newsinstereo = 0;
6487 newerror[0] = std::sqrt(
covmat(0, 0));
6489 if (state_maxsipull->
sinStereo() != 0) {
6507 newerror[0] = std::sqrt(
v0);
6510 double cosstereo = (newsinstereo == 0) ? 1. : std::sqrt(1 - newsinstereo * newsinstereo);
6512 if (cosstereo != 1.) {
6514 cosstereo * (broadrot->localParameters()[
Trk::locX] - trackpar_maxsipull->parameters()[
Trk::locX]) +
6515 newsinstereo * (broadrot->localParameters()[
Trk::locY] - trackpar_maxsipull->parameters()[
Trk::locY])
6518 newres1 = broadrot->localParameters()[
Trk::locX] - trackpar_maxsipull->parameters()[
Trk::locX];
6521 if (newerror[0] == 0.0) {
6522 ATH_MSG_WARNING(
"Measurement error is zero or negative, treating as outlier");
6526 newpull1 = std::abs(newres1 / newerror[0]);
6530 newerror[1] = std::sqrt(
covmat(1, 1));
6531 newres2 = broadrot->localParameters()[
Trk::locY] - trackpar_maxsipull->parameters()[
Trk::locY];
6532 newpull2 = std::abs(newres2 / newerror[1]);
6535 newpull =
std::max(newpull1, newpull2);
6541 (newerror[0] > 1.5 * olderror[0] || newerror[1] > 1.5 * std::abs(olderror[1]))
6543 if ((measno_maxsipull < 0) or(measno_maxsipull >= (
int)
res.size())) {
6544 throw std::runtime_error(
6545 "'res' array index out of range in TrkGlobalChi2Fitter/src/GlobalChi2Fitter.cxx:" +
std::to_string(__LINE__)
6550 newtrajectory = oldtrajectory;
6552 if (
a.cols() != nfitpars) {
6556 double oldres1 =
res[measno_maxsipull];
6557 res[measno_maxsipull] = newres1;
6558 err[measno_maxsipull] = newerror[0];
6560 for (
int i = 0;
i < nfitpars;
i++) {
6561 if (weightderiv(measno_maxsipull,
i) == 0) {
6565 b[
i] -= weightderiv(measno_maxsipull,
i) * (oldres1 / olderror[0] - (newres1 * olderror[0]) / (newerror[0] * newerror[0]));
6567 for (
int j =
i; j < nfitpars; j++) {
6571 weightderiv(measno_maxsipull,
i) *
6572 weightderiv(measno_maxsipull, j) *
6573 ((olderror[0] * olderror[0]) / (newerror[0] * newerror[0]) - 1)
6577 weightderiv(measno_maxsipull,
i) *= olderror[0] / newerror[0];
6581 double oldres2 =
res[measno_maxsipull + 1];
6582 res[measno_maxsipull + 1] = newres2;
6583 err[measno_maxsipull + 1] = newerror[1];
6585 for (
int i = 0;
i < nfitpars;
i++) {
6586 if (weightderiv(measno_maxsipull + 1,
i) == 0) {
6590 b[
i] -= weightderiv(measno_maxsipull + 1,
i) * (oldres2 / olderror[1] - (newres2 * olderror[1]) / (newerror[1] * newerror[1]));
6592 for (
int j =
i; j < nfitpars; j++) {
6596 weightderiv(measno_maxsipull + 1,
i) *
6597 weightderiv(measno_maxsipull + 1, j) *
6598 ((olderror[1] * olderror[1]) / (newerror[1] * newerror[1]) - 1)
6603 weightderiv(measno_maxsipull + 1,
i) *= olderror[1] / newerror[1];
6608 "Recovering outlier, hitno=" << hitno_maxsipull <<
" measno=" <<
6609 measno_maxsipull <<
" pull=" << maxsipull <<
" olderror_0=" <<
6610 olderror[0] <<
" newerror_0=" << newerror[0] <<
" olderror_1=" <<
6611 olderror[1] <<
" newerror_1=" << newerror[1]
6620 ((n3sigma < 2 && maxsipull > cut2 && maxsipull <
cut) || n3sigma > 1) &&
6630 "Removing outlier, hitno=" << hitno_maxsipull <<
", measno=" <<
6631 measno_maxsipull <<
" pull=" << maxsipull
6638 cleanup_newtrajectory = std::make_unique<GXFTrajectory>(*oldtrajectory);
6639 newtrajectory = cleanup_newtrajectory.get();
6641 if (newa.cols() != nfitpars) {
6647 if ((measno_maxsipull < 0) or(measno_maxsipull >= (
int)
res.size())) {
6648 throw std::runtime_error(
6649 "'res' array index out of range in TrkGlobalChi2Fitter/src/GlobalChi2Fitter.cxx:" +
std::to_string(__LINE__)
6653 double oldres1 =
res[measno_maxsipull];
6654 newres[measno_maxsipull] = 0;
6656 for (
int i = 0;
i < nfitpars;
i++) {
6657 if (weightderiv(measno_maxsipull,
i) == 0) {
6661 newb[
i] -= weightderiv(measno_maxsipull,
i) * oldres1 / olderror[0];
6663 for (
int j =
i; j < nfitpars; j++) {
6667 weightderiv(measno_maxsipull,
i) *
6668 weightderiv(measno_maxsipull, j)
6672 newweightderiv(measno_maxsipull,
i) = 0;
6676 double oldres2 =
res[measno_maxsipull + 1];
6677 newres[measno_maxsipull + 1] = 0;
6679 for (
int i = 0;
i < nfitpars;
i++) {
6680 if (weightderiv(measno_maxsipull + 1,
i) == 0) {
6684 newb[
i] -= weightderiv(measno_maxsipull + 1,
i) * oldres2 / olderror[1];
6686 for (
int j =
i; j < nfitpars; j++) {
6687 if (weightderiv(measno_maxsipull + 1, j) == 0) {
6694 weightderiv(measno_maxsipull + 1,
i) *
6695 weightderiv(measno_maxsipull + 1, j)
6699 newweightderiv(measno_maxsipull + 1,
i) = 0;
6703 newtrajectory->
setOutlier(stateno_maxsipull);
6727 ctx, cache, *newtrajectory,
it, *newap, *newbp, lu_m, doderiv);
6743 double oldchi2 = oldtrajectory->
chi2() / oldtrajectory->
nDOF();
6744 double newchi2 = (newtrajectory->
nDOF() > 0) ? newtrajectory->
chi2() / newtrajectory->
nDOF() : 0;
6747 if (newtrajectory->
nDOF() != oldtrajectory->
nDOF() && maxsipull > cut2) {
6748 mindiff = (oldchi2 > .33 *
m_chi2cut || noutl > 0) ? .8 : 1.;
6750 if (noutl == 0 && maxsipull <
cut - .5 && oldchi2 < .5 *
m_chi2cut) {
6755 if (newchi2 > oldchi2 || (newchi2 > oldchi2 - mindiff && newchi2 > .33 * oldchi2)) {
6756 ATH_MSG_DEBUG(
"Outlier not confirmed, keeping old trajectory");
6764 (void)cleanup_oldtrajectory.release();
6765 return oldtrajectory;
6767 if (oldtrajectory != newtrajectory) {
6768 cleanup_oldtrajectory = std::move(cleanup_newtrajectory);
6769 oldtrajectory = newtrajectory;
6776 Eigen::LLT < Eigen::MatrixXd > lltOfW(
a);
6777 if (lltOfW.info() == Eigen::Success) {
6781 int ncols =
a.cols();
6782 Amg::MatrixX weightInvAMG = Amg::MatrixX::Identity(ncols, ncols);
6783 fullcov = lltOfW.solve(weightInvAMG);
6801 oldtrajectory->
nDOF() > 0 &&
6810 (void)cleanup_oldtrajectory.release();
6811 return oldtrajectory;
6823 if (
const auto *pMeas{hit->measurement()};
6829 nrealmeas += hit->numberOfMeasuredParameters();
6839 if (
const auto *pMeas{hit->measurement()};
6845 for (
int i = measindex;
i < measindex + hit->numberOfMeasuredParameters();
i++) {
6848 if ((j == 4 && !oldtrajectory.
m_straightline) || j >= nperpars + 2 * nscat) {
6856 measindex += hit->numberOfMeasuredParameters();
6857 }
else if (hit->materialEffects() ==
nullptr) {
6858 measindex2 += hit->numberOfMeasuredParameters();
6864 const EventContext & ctx,
6869 GXFTrackState *firstmeasstate =
nullptr, *lastmeasstate =
nullptr;
6871 std::unique_ptr<const TrackParameters> per(
nullptr);
6874 std::unique_ptr<const TrackParameters> prevpar(
6879 std::vector<std::pair<const Layer *, const Layer *>> & upstreamlayers = oldtrajectory.
upstreamMaterialLayers();
6882 for (
int i = (
int)upstreamlayers.size() - 1;
i >= 0;
i--) {
6883 if (prevpar ==
nullptr) {
6890 if (
layer ==
nullptr) {
6891 layer = upstreamlayers[
i].second;
6895 prevpar->position(), prevpar->momentum().unit()
6913 std::unique_ptr<const TrackParameters> layerpar(
6917 layer->surfaceRepresentation(),
6925 if (layerpar ==
nullptr) {
6929 if (
layer->surfaceRepresentation().bounds().inside(layerpar->localPosition())) {
6933 prevpar = std::move(layerpar);
6946 if (startfactor > 0.5) {
6947 std::unique_ptr<const TrackParameters> updatedpar =
m_matupdator->update(
6951 if (updatedpar !=
nullptr) {
6961 const Layer *endlayer = lastmeasstate->trackParameters()->associatedSurface().associatedLayer();
6971 if (endfactor > 0.5) {
6972 std::unique_ptr<const TrackParameters> updatedpar =
m_matupdator->update(
6973 lastmeasstate->trackParameters(), *endlayer,
alongMomentum, matEffects
6976 if (updatedpar !=
nullptr) {
6977 lastmeasstate->setTrackParameters(std::move(updatedpar));
6982 if (prevpar !=
nullptr) {
6994 if (per ==
nullptr) {
6995 ATH_MSG_DEBUG(
"Failed to extrapolate to perigee, returning 0");
7014 std::unique_ptr<GXFTrackState>
7016 const EventContext & ctx,
7023 if (per ==
nullptr) {
7027 ATH_MSG_DEBUG(
"Final perigee: " << *per <<
" pos: " << per->position() <<
" pT: " << per->pT());
7033 const std::vector<std::unique_ptr<TrackParameters>> & hc,
7034 std::set<Identifier> & id_set,
7035 std::set<Identifier> & sct_set,
7045 for (
const std::unique_ptr<TrackParameters> &
tp : hc) {
7051 if (
tp ==
nullptr) {
7062 if (de ==
nullptr) {
7073 if (id_set.find(
id) != id_set.end()) {
7122 if (sct_set.find(
os) != sct_set.end()) {
7123 ++
rv.m_sct_double_hole;
7152 for (
const std::unique_ptr<GXFTrackState> &
s : trajectory.
trackStates()) {
7164 std::vector<std::reference_wrapper<GXFTrackState>>
rv;
7172 for (
const std::unique_ptr<GXFTrackState> &
s : trajectory.
trackStates()) {
7186 rv.emplace_back(*
s);
7193 const TrkDetElementBase * de =
s->trackParameters()->associatedSurface().associatedDetectorElement();
7195 if (de !=
nullptr) {
7208 if (
s.get() == lastmeas) {
7218 const EventContext & ctx,
7219 const std::vector<std::reference_wrapper<GXFTrackState>> &
states
7231 constexpr
uint min_meas = 3;
7232 if (std::count_if(
states.begin(),
states.end(), [](
const GXFTrackState &
s){ return s.getStateType(TrackStateOnSurface::Measurement); }) < min_meas) {
7236 bool seen_meas =
false;
7238 std::set<Identifier> id_set;
7239 std::set<Identifier> sct_set;
7245 for (std::size_t
i = 0;
i <
states.size() - 1;
i++) {
7268 double dist = (
beg.trackParameters()->position() -
end.trackParameters()->position()).
norm();
7270 bool zStartValid = std::abs(
beg.trackParameters()->position().z())<10000.;
7272 ATH_MSG_DEBUG(
"Pathological track parameter well outside of detector");
7273 ATH_MSG_DEBUG(
"Propagator might have issue with this, skipping");
7283 if (seen_meas && dist >= 2.5 && zStartValid) {
7290 std::optional<std::vector<std::unique_ptr<TrackParameters>>> & hc =
beg.getHoles();
7291 std::vector<std::unique_ptr<TrackParameters>>
states;
7299 if (hc.has_value()) {
7338 std::vector<std::unique_ptr<Trk::TrackParameters>> bl =
m_extrapolator->extrapolateBlindly(
7359 const EventContext & ctx,
7365 auto trajectory = std::make_unique<Trk::TrackStates>();
7373 std::unique_ptr<GXFTrackState> perigee_ts =
makeTrackFindPerigee(ctx, cache, oldtrajectory, matEffects);
7375 if (perigee_ts ==
nullptr) {
7381 trajectory->reserve(tmptrajectory.
trackStates().size());
7387 hit->resetTrackCovariance();
7393 hit->materialEffects()))
7399 auto trackState = hit->trackStateOnSurface();
7400 hit->resetTrackCovariance();
7401 trajectory->emplace_back(trackState.release());
7404 auto qual = std::make_unique<FitQuality>(tmptrajectory.
chi2(), tmptrajectory.
nDOF());
7424 std::unique_ptr<Track>
rv = std::make_unique<Track>(
info, std::move(trajectory), std::move(
qual));
7434 std::unique_ptr<TrackSummary>
ts = std::make_unique<TrackSummary>();
7443 std::optional<TrackHoleCount> hole_count;
7468 if (hole_count.has_value()) {
7481 rv->setTrackSummary(std::move(
ts));
7487 GlobalChi2Fitter::~GlobalChi2Fitter() =
default;
7490 const EventContext & ctx,
7499 std::vector<std::unique_ptr<TrackParameters>>
rv =
m_extrapolator->extrapolateStepwise(
7513 &
rv.front()->associatedSurface() == &
src.associatedSurface() ||
7514 trackParametersClose(*
rv.front(),
src, 0.001) ||
7518 rv.front().reset(
nullptr);
7529 &
rv.back()->associatedSurface() == &
src.associatedSurface() ||
7530 trackParametersClose(*
rv.back(),
src, 0.001) ||
7534 rv.back().reset(
nullptr);
7541 const EventContext & ctx,
7549 std::unique_ptr<const TrackParameters>
rv;
7550 std::optional<TransportJacobian> jac{};
7561 if (
rv !=
nullptr && calcderiv) {
7566 std::optional<std::vector<std::unique_ptr<TrackParameters>>> extrapolation;
7575 std::move(extrapolation)
7580 const EventContext & ctx,
7591 ctx, prev,
ts, propdir, bf, calcderiv, holesearch
7594 if (
rv.m_parameters ==
nullptr) {
7595 propdir = invertPropdir(propdir);
7598 ctx, prev,
ts, propdir, bf, calcderiv, holesearch
7606 const EventContext& ctx,
7613 std::vector<std::unique_ptr<GXFTrackState>> &
states = trajectory.
trackStates();
7616 std::unique_ptr<const TrackParameters> tmptrackpar;
7618 for (
int hitno = nstatesupstream - 1; hitno >= 0; hitno--) {
7648 (
rv.m_parameters !=
nullptr) &&
7649 (prevtrackpar->
position() -
rv.m_parameters->position()).mag() > 5 *
mm
7655 if (
rv.m_parameters ==
nullptr) {
7656 ATH_MSG_DEBUG(
"propagation failed, prev par: " << *prevtrackpar <<
7657 " pos: " << prevtrackpar->
position() <<
" destination surface: " << surf1);
7661 states[hitno]->setTrackParameters(std::move(
rv.m_parameters));
7665 if (
rv.m_jacobian != std::nullopt) {
7667 states[hitno]->materialEffects() !=
nullptr &&
7668 states[hitno]->materialEffects()->deltaE() != 0 &&
7669 states[hitno]->materialEffects()->sigmaDeltaE() <= 0 &&
7672 double p = 1. / std::abs(currenttrackpar->parameters()[
Trk::qOverP]);
7673 double de = std::abs(
states[hitno]->materialEffects()->deltaE());
7675 double newp = std::sqrt(
p *
p + 2 * de * std::sqrt(
mass *
mass +
p *
p) + de * de);
7676 (*
rv.m_jacobian) (4, 4) = ((
p +
p * de / std::sqrt(
p *
p +
mass *
mass)) / newp) *
p *
p / (newp * newp);
7679 states[hitno]->setJacobian(*
rv.m_jacobian);
7680 }
else if (calcderiv) {
7687 if (meff !=
nullptr && hitno != 0) {
7689 surf, *meff, *
states[hitno]->trackParameters(), trajectory.
mass(), -1
7692 if (std::holds_alternative<FitterStatusCode>(
r)) {
7693 return std::get<FitterStatusCode>(
r);
7696 tmptrackpar = std::move(
std::get<std::unique_ptr<const TrackParameters>>(
r));
7697 prevtrackpar = tmptrackpar.get();
7699 prevtrackpar = currenttrackpar;
7705 for (
int hitno = nstatesupstream; hitno < (
int)
states.size(); hitno++) {
7727 (
rv.m_parameters !=
nullptr) &&
7729 (prevtrackpar->
position() -
rv.m_parameters->position()).mag() > 5 *
mm
7734 if (
rv.m_parameters ==
nullptr) {
7735 ATH_MSG_DEBUG(
"propagation failed, prev par: " << *prevtrackpar <<
7736 " pos: " << prevtrackpar->
7737 position() <<
" destination surface: " << surf);
7741 if (
rv.m_jacobian != std::nullopt) {
7743 states[hitno]->materialEffects() !=
nullptr &&
7744 states[hitno]->materialEffects()->deltaE() != 0 &&
7745 states[hitno]->materialEffects()->sigmaDeltaE() <= 0 &&
7748 double p = 1 / std::abs(
rv.m_parameters->parameters()[
Trk::qOverP]);
7749 double de = std::abs(
states[hitno]->materialEffects()->deltaE());
7751 double newp =
p *
p - 2 * de * std::sqrt(
mass *
mass +
p *
p) + de * de;
7754 newp = std::sqrt(newp);
7757 (*
rv.m_jacobian) (4, 4) = ((
p -
p * de / std::sqrt(
p *
p +
mass *
mass)) / newp) *
p *
p / (newp * newp);
7760 states[hitno]->setJacobian(*
rv.m_jacobian);
7761 }
else if (calcderiv) {
7768 if (meff !=
nullptr) {
7770 surf, *meff, *
rv.m_parameters, trajectory.
mass(), +1
7773 if (std::holds_alternative<FitterStatusCode>(
r)) {
7774 return std::get<FitterStatusCode>(
r);
7777 rv.m_parameters = std::move(
std::get<std::unique_ptr<const TrackParameters>>(
r));
7780 states[hitno]->setTrackParameters(std::move(
rv.m_parameters));
7781 prevtrackpar =
states[hitno]->trackParameters();
7800 ATH_MSG_DEBUG(
"Angles out of range, phi: " << newphi <<
" theta: " << newtheta);
7804 double newqoverp = 0;
7818 newqoverp = std::copysign(1 / std::sqrt(newp2),
old[
Trk::qOverP]);
7825 old[0],
old[1], newphi, newtheta, newqoverp, std::nullopt
7837 using Matrix55 = Eigen::Matrix<double, 5, 5>;
7839 Matrix55 initialjac;
7840 initialjac.setZero();
7841 initialjac(4, 4) = 1;
7843 Matrix55 jacvertex(initialjac);
7845 std::vector<Matrix55, Eigen::aligned_allocator<Matrix55>> jacscat(trajectory.
numberOfScatterers(), initialjac);
7846 std::vector<Matrix55, Eigen::aligned_allocator<Matrix55>> jacbrem(trajectory.
numberOfBrems(), initialjac);
7848 std::vector<std::unique_ptr<GXFTrackState>> &
states = trajectory.
trackStates();
7851 int hit_begin = 0, hit_end = 0, scatno = 0, bremno = 0;
7853 for (
bool forward : {
false,
true}) {
7855 hit_begin = nstatesupstream;
7857 scatno = nscatupstream;
7858 bremno = nbremupstream;
7860 hit_begin = nstatesupstream - 1;
7867 int hitno = hit_begin;
7868 forward ? (hitno < hit_end) : (hitno >= hit_end);
7869 hitno += (forward ? 1 : -1)
7872 state =
states[hitno].get();
7876 if (fillderivmat && state->derivatives().cols() != nfitpars) {
7877 state->derivatives().resize(5, nfitpars);
7878 state->derivatives().setZero();
7881 int jminscat = 0, jmaxscat = 4, jminbrem = 0, jmaxbrem = 4;
7883 if (hitno == (forward ? hit_end - 1 : 0)) {
7884 if (!fillderivmat) {
7892 Eigen::Matrix<double, 5, 5> & jac = state->jacobian();
7894 if (hitno == nstatesupstream + (forward ? 0 : -1)) {
7895 jacvertex.block<4, 5>(0, 0) = jac.block<4, 5>(0, 0);
7896 jacvertex(4, 4) = jac(4, 4);
7898 int jmin = 0, jmax = 0, jcnt = 0;
7899 int lp_bgn = 0, lp_end = 0;
7903 jcnt = jmax - jmin + 1;
7905 lp_bgn = forward ? nscatupstream : nscatupstream - 1;
7908 for (
int i = lp_bgn; forward ? (
i < lp_end) : (
i > lp_end);
i += (forward ? 1 : -1)) {
7910 i == scatno + (forward ? -1 : 1) &&
7911 prevstate !=
nullptr &&
7915 jacscat[
i].block(0, jmin, 4, jcnt) = jac.block(0, jmin, 4, jcnt);
7916 jacscat[
i](4, 4) = jac(4, 4);
7918 calculateJac(jac, jacscat[
i], jmin, jmax);
7922 Eigen::MatrixXd & derivmat = state->derivatives();
7923 int scatterPos = nperpars + 2 *
i;
7925 derivmat.block<4, 2>(0, scatterPos) = (forward ? 1 : -1) * jacscat[
i].block<4, 2>(0, 2);
7931 jcnt = jmax - jmin + 1;
7933 lp_bgn = forward ? nbremupstream : nbremupstream - 1;
7936 for (
int i = lp_bgn; forward ? (
i < lp_end) : (
i > lp_end);
i += (forward ? 1 : -1)) {
7938 i == bremno + (forward ? -1 : 1) &&
7943 jacbrem[
i].block(0, jmin, 4, jcnt) = jac.block(0, jmin, 4, jcnt);
7944 jacbrem[
i](4, 4) = jac(4, 4);
7946 calculateJac(jac, jacbrem[
i], jmin, jmax);
7950 Eigen::MatrixXd & derivmat = state->derivatives();
7951 int scatterPos = nperpars + 2 * nscats +
i;
7953 derivmat.block<5, 1>(0, scatterPos) = (forward ? .001 : -.001) * jacbrem[
i].block<5, 1>(0, 4);
7957 calculateJac(jac, jacvertex, 0, 4);
7961 Eigen::MatrixXd & derivmat = state->derivatives();
7962 derivmat.block(0, 0, 4, nperpars) = jacvertex.block(0, 0, 4, nperpars);
7964 if (nperpars == 5) {
7965 derivmat.col(4).segment(0, 4) *= .001;
7966 derivmat(4, 4) = .001 * jacvertex(4, 4);
7972 (!trajectory.
prefit() ||
states[hitno]->materialEffects()->deltaE() == 0)
7974 scatno += (forward ? 1 : -1);
7978 states[hitno]->materialEffects() &&
7979 states[hitno]->materialEffects()->sigmaDeltaE() > 0
7981 bremno += (forward ? 1 : -1);
7984 prevstate =
states[hitno].get();
7993 bool onlylocal)
const {
7999 std::vector<std::unique_ptr<GXFTrackState>> &
states = trajectory.
trackStates();
8003 int i = nstatesupstream;
8004 for (
int j = 0; j < (
int)
states.size(); j++) {
8005 if (j < nstatesupstream) {
8012 for (
int stateno = 0; stateno < (
int)
states.size(); stateno++) {
8013 if (stateno == 0 || stateno == nstatesupstream) {
8014 prevstate =
nullptr;
8017 std::unique_ptr<GXFTrackState> & state =
states[
index];
8018 if (state->materialEffects() !=
nullptr) {
8019 prevstate = state.get();
8023 if (!state->hasTrackCovariance()) {
8024 state->zeroTrackCovariance();
8026 AmgMatrix(5, 5) & trackerrmat = state->trackCovariance();
8028 if ((prevstate !=
nullptr) &&
8032 Eigen::Matrix<double, 5, 5> & jac = state->jacobian();
8035 trackerrmat = jac * prevcov * jac.transpose();
8039 trackerrmat = derivatives * fullcovmat * derivatives.transpose();
8050 bool errorok =
true;
8051 for (
int i = 0;
i < 5;
i++) {
8054 && trackerrmat(
i,
i) > meascov(j, j)) {
8056 double scale = std::sqrt(meascov(j, j) / trackerrmat(
i,
i));
8057 trackerrmat(
i,
i) = meascov(j, j);
8058 for (
int k = 0;
k < 5;
k++) {
8068 for (
int i = 0;
i < 5;
i++) {
8072 for (
int j = 0; j < 5; j++) {
8080 trackerrmat(4, 4) = 1
e-20;
8084 state->trackParameters();
8086 std::optional<
AmgMatrix(5, 5)> trkerrmat;
8088 if (state->hasTrackCovariance()) {
8089 trkerrmat = (state->trackCovariance());
8091 trkerrmat = std::nullopt;
8094 const AmgVector(5) & tpars = tmptrackpar->parameters();
8095 std::unique_ptr<const TrackParameters> trackpar(
8101 std::move(trkerrmat))
8103 state->setTrackParameters(std::move(trackpar));
8106 if (errorok && trajectory.
nDOF() > 0) {
8107 fitQual =
m_updator->fullStateFitQuality(
8108 *state->trackParameters(),
8116 state->setFitQuality(fitQual);
8118 prevstate = state.get();
8122 std::optional<TransportJacobian>
8124 const EventContext& ctx,
8138 std::optional<TransportJacobian> jac = std::make_optional<TransportJacobian>(J);
8141 0.01, 0.01, 0.00001, 0.00001, 0.000000001
8152 for (
int i = 0;
i < 5;
i++) {
8155 if (thisdiscsurf &&
i == 1) {
8159 vecpluseps[paraccessor.
pardef[
i]] += eps[
i];
8160 vecminuseps[paraccessor.
pardef[
i]] -= eps[
i];
8161 if (thiscylsurf &&
i == 0) {
8162 if (vecpluseps[0] / previousSurface.
bounds().
r() >
M_PI) {
8163 vecpluseps[0] -= 2 *
M_PI * previousSurface.
bounds().
r();
8165 if (vecminuseps[0] / previousSurface.
bounds().
r() < -
M_PI) {
8166 vecminuseps[0] += 2 *
M_PI * previousSurface.
bounds().
r();
8169 if (thisdiscsurf &&
i == 1) {
8170 if (vecpluseps[
i] >
M_PI) {
8171 vecpluseps[
i] -= 2 *
M_PI;
8173 if (vecminuseps[
i] < -
M_PI) {
8174 vecminuseps[
i] += 2 *
M_PI;
8180 std::unique_ptr<const TrackParameters> parpluseps(
8190 std::unique_ptr<const TrackParameters> parminuseps(
8201 std::unique_ptr<const TrackParameters> newparpluseps(
8212 std::unique_ptr<const TrackParameters> newparminuseps(
8227 if (newparpluseps ==
nullptr) {
8239 if (newparminuseps ==
nullptr) {
8251 if ((newparpluseps ==
nullptr) || (newparminuseps ==
nullptr)) {
8255 for (
int j = 0; j < 5; j++) {
8256 double diff = newparpluseps->parameters()[paraccessor.
pardef[j]] -
8257 newparminuseps->parameters()[paraccessor.
pardef[j]];
8258 if (cylsurf && j == 0) {
8268 if (discsurf && j == 1) {
8269 if (std::abs(std::abs(
diff) - 2 *
M_PI) < std::abs(
diff)) {
8278 (*jac) (j,
i) =
diff / (2 * eps[
i]);
8291 (
"Configure the minimum number of Iterations via jobOptions");
8316 auto nmeas1 = pDataVector->size();
8317 const auto *pLastValue = (*pDataVector)[nmeas1 - 1];
8325 testrot =
static_cast<const RIO_OnTrack *
>(pLastValue);
8327 if (lastMeasIsCompetingRIO){
8333 if (testrot ==
nullptr) {
8334 const auto *pPenultimate = (*pDataVector)[nmeas1 - 2];
8337 if(penultimateIsRIO){
8338 testrot =
static_cast<const RIO_OnTrack *
>(pPenultimate);
8340 if (penultimateIsCompetingRIO){
8349 (testrot !=
nullptr) &&
8366 if (cond_obj ==
nullptr) {
8375 std::stringstream
msg;
8377 throw std::runtime_error(
msg.str());