56 #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 const double limited = std::clamp(magnitude, lo, hi);
103 return std::copysign(limited,
qOverP);
107 std::pair<const Trk::TrackParameters *, const Trk::TrackParameters *> getFirstLastIdPar(
const Trk::Track &
track) {
113 while ((firstidpar ==
nullptr) && parit !=
track.trackParameters()->end()) {
115 ((**parit).covariance() !=
nullptr) &&
124 parit =
track.trackParameters()->end();
128 ((**parit).covariance() !=
nullptr) &&
133 }
while ((lastidpar ==
nullptr) && parit !=
track.trackParameters()->begin());
135 return std::make_pair(firstidpar, lastidpar);
150 std::abs(
a.parameters()[0] -
b.parameters()[0]) <
e &&
151 std::abs(
a.parameters()[1] -
b.parameters()[1]) <
e &&
152 std::abs(
a.parameters()[2] -
b.parameters()[2]) <
e
163 calculateJac(Eigen::Matrix<double, 5, 5> &jac,
164 Eigen::Matrix<double, 5, 5> &
out,
165 int jmin,
int jmax) {
169 out.block(0, 0, 4, jmin).setZero();
173 out.block(0, jmax + 1, 4, 5 - (jmax + 1)).setZero();
176 out(4, 4) = jac(4, 4);
186 std::pair<double, double> principalComponentAnalysis2x2(
const Amg::MatrixX &
mat) {
187 const double trace =
mat(0, 0) +
mat(1, 1);
188 const double diagonalProduct =
mat(0, 0) *
mat(1, 1);
189 const double mat01Sq =
mat(0, 1) *
mat(0, 1);
190 const double discriminant = std::sqrt(trace * trace - 4. * (diagonalProduct - mat01Sq));
192 const double eigenValueSmall = 0.5 * (trace -
discriminant);
193 const double stereoAngle = 0.5 * std::asin(2 *
mat(0, 1) / (-
discriminant));
195 return std::make_pair(eigenValueSmall, stereoAngle);
201 const std::string &
t,
202 const std::string &
n,
230 ATH_MSG_ERROR(
"Hole search requested but no boundary check tool provided.");
231 return StatusCode::FAILURE;
256 ATH_MSG_WARNING(
"FillDerivativeMatrix option selected, switching off acceleration!");
280 ATH_MSG_ERROR(
"Hole search requested but track summaries are disabled.");
281 return StatusCode::FAILURE;
286 return StatusCode::SUCCESS;
292 if (m_fit_status[
S_FITS] > 0) {
295 <<
" track fits failed because of a matrix inversion failure");
297 <<
" tracks were rejected by the outlier logic");
299 <<
" track fits failed because of a propagation failure");
301 <<
" track fits failed because of an invalid angle (theta/phi)");
303 <<
" track fits failed because the fit did not converge");
305 <<
" tracks did not pass the chi^2 cut");
307 <<
" tracks were killed by the energy loss update");
310 return StatusCode::SUCCESS;
315 std::unique_ptr<Track>
317 const EventContext& ctx,
323 ATH_MSG_DEBUG(
"--> entering GlobalChi2Fitter::fit(Track,Track,)");
338 const Track *indettrack = firstismuon ? &intrk2 : &intrk1;
339 const Track *muontrack = firstismuon ? &intrk1 : &intrk2;
341 bool measphi =
false;
348 rot = &crot->rioOnTrack(0);
358 const double dotprod2 = measdir.dot(
361 if (std::abs(dotprod1) < .5 && std::abs(dotprod2) < .5) {
369 auto [firstidpar, lastidpar] = getFirstLastIdPar(*indettrack);
371 if ((firstidpar ==
nullptr) || (lastidpar ==
nullptr)) {
375 std::unique_ptr<const TrackParameters> parforcalo =
unique_clone(firstismuon ? firstidpar : lastidpar);
378 const AmgVector(5) & newpars = parforcalo->parameters();
380 parforcalo=parforcalo->associatedSurface().createUniqueTrackParameters(
381 newpars[0], newpars[1], newpars[2], newpars[3], 1 / 5000., std::nullopt
385 std::vector < MaterialEffectsOnTrack > calomeots;
388 calomeots =
m_calotool->extrapolationSurfacesAndEffects(
392 parforcalo->associatedSurface(),
405 if (calomeots.empty()) {
410 std::unique_ptr<Track>
track;
415 const bool tmp4 = cache.
m_idmat;
420 const double qoverpid = measperid !=
nullptr ? measperid->parameters()[
Trk::qOverP] : 0;
421 const double qoverpmuon = measpermuon !=
nullptr ? measpermuon->parameters()[
Trk::qOverP] : 0;
423 const AmgSymMatrix(5) * errmatid = measperid !=
nullptr ? measperid->covariance() :
nullptr;
424 const AmgSymMatrix(5) * errmatmuon = measpermuon !=
nullptr ? measpermuon->covariance() :
nullptr;
428 (errmatid !=
nullptr) &&
429 (errmatmuon !=
nullptr) &&
434 const double piderror = std::sqrt((*errmatid) (4, 4)) / (qoverpid * qoverpid);
435 const double pmuonerror = std::sqrt((*errmatmuon) (4, 4)) / (qoverpmuon * qoverpmuon);
436 const double energyerror = std::sqrt(
437 calomeots[1].energyLoss()->sigmaDeltaE() *
438 calomeots[1].energyLoss()->sigmaDeltaE() + piderror * piderror +
439 pmuonerror * pmuonerror
443 (std::abs(calomeots[1].energyLoss()->deltaE()) -
444 std::abs(1 / qoverpid) + std::abs(1 / qoverpmuon)
447 ATH_MSG_DEBUG(
"Changing from measured to parametrized energy loss");
452 parforcalo->associatedSurface(),
457 if (calomeots.empty()) {
465 bool firstfitwasattempted =
false;
468 if (!caloEntranceIsValid) {
479 qoverpid * qoverpmuon > 0
485 firstfitwasattempted =
true;
490 (
track ==
nullptr) &&
491 !firstfitwasattempted &&
498 trajectory = trajectory2;
502 bool pseudoupdated =
false;
504 if (
track !=
nullptr) {
505 for (std::unique_ptr<GXFTrackState> & pseudostate : trajectory.
trackStates()) {
506 if (pseudostate ==
nullptr) {
517 if ((pseudostate ==
nullptr) || pseudostate->fitQuality().chiSquared() < 10) {
522 const std::unique_ptr<const TrackParameters> updpar(
m_updator->removeFromState(
524 pseudostate->measurement()->localParameters(),
525 pseudostate->measurement()->localCovariance()
528 if (updpar ==
nullptr) {
535 std::unique_ptr<const PseudoMeasurementOnTrack> newpseudo = std::make_unique<const PseudoMeasurementOnTrack>(
543 pseudostate->setMeasurement(std::move(newpseudo));
547 pseudostate->setMeasurementErrors(
errors);
548 pseudoupdated =
true;
559 *
track->perigeeParameters(),
570 if (
track !=
nullptr) {
571 track->info().addPatternReco(intrk1.
info());
572 track->info().addPatternReco(intrk2.
info());
583 const EventContext& ctx,
585 const Track & intrk1,
586 const Track & intrk2,
588 std::vector<MaterialEffectsOnTrack> & calomeots
590 ATH_MSG_DEBUG(
"--> entering GlobalChi2Fitter::mainCombinationStrategy");
595 const Track *indettrack = firstismuon ? &intrk2 : &intrk1;
596 const Track *muontrack = firstismuon ? &intrk1 : &intrk2;
598 auto [tmpfirstidpar, tmplastidpar] = getFirstLastIdPar(*indettrack);
599 std::unique_ptr<const TrackParameters> firstidpar =
unique_clone(tmpfirstidpar);
600 std::unique_ptr<const TrackParameters> lastidpar =
unique_clone(tmplastidpar);
602 if ((firstidpar ==
nullptr) || (lastidpar ==
nullptr)) {
616 std::unique_ptr<const TrackParameters> tp_closestmuon =
nullptr;
618 while (closestmuonmeas ==
nullptr) {
619 closestmuonmeas =
nullptr;
622 if ((**tsosit).measurementOnTrack() !=
nullptr) {
623 closestmuonmeas = (**tsosit).measurementOnTrack();
625 if (thispar !=
nullptr) {
626 const AmgVector(5) & parvec = thispar->parameters();
628 parvec[0], parvec[1], parvec[2], parvec[3], parvec[4], std::nullopt
642 std::unique_ptr<const TrackParameters> tmppar;
645 if ((tp_closestmuon !=
nullptr) && msEntranceIsValid) {
650 std::unique_ptr<const std::vector<const TrackStateOnSurface *>> matvec;
652 if (tmppar !=
nullptr) {
653 const Surface & associatedSurface = tmppar->associatedSurface();
654 std::unique_ptr<Surface> muonsurf =
nullptr;
660 const double radius = cylbounds->
r();
662 muonsurf = std::make_unique<CylinderSurface>(trans,
radius + 1, hlength);
666 const double newz = (
667 associatedSurface.
center().z() > 0 ?
668 associatedSurface.
center().z() + 1 :
669 associatedSurface.
center().z() - 1
673 associatedSurface.
center().x(),
674 associatedSurface.
center().y(),
678 trans.translation() << newpos;
681 const double rmin = discbounds->
rMin();
682 const double rmax = discbounds->
rMax();
683 muonsurf = std::make_unique<DiscSurface>(trans, rmin, rmax);
687 if (muonsurf !=
nullptr) {
699 std::vector<const TrackStateOnSurface *> tmp_matvec;
701 if ((matvec !=
nullptr) && !matvec->empty()) {
702 tmp_matvec = *matvec;
703 delete tmp_matvec.back();
704 tmp_matvec.pop_back();
706 for (
auto &
i : tmp_matvec) {
725 if (tmppar ==
nullptr) {
739 if (tmppar ==
nullptr) {
743 AmgVector(5) newpars = tmppar->parameters();
748 const double oldp = std::abs(1 / newpars[
Trk::qOverP]);
749 const double newp2 = oldp * oldp + (!firstismuon ? 2 : -2) * de * std::sqrt(
mass *
mass + oldp * oldp) + de * de;
756 tp_closestmuon=tmppar->associatedSurface().createUniqueTrackParameters(
757 newpars[0], newpars[1], newpars[2], newpars[3], newpars[4], std::nullopt
773 for (; itStates != endState; ++itStates) {
780 if ((*itStates)->materialEffectsOnTrack() !=
nullptr) {
787 const auto *
const pBaseMEOT = (*itStates)->materialEffectsOnTrack();
791 const auto *
const pMEOT =
static_cast<const MaterialEffectsOnTrack *
>((*itStates)->materialEffectsOnTrack());
792 if ((pMEOT->scatteringAngles() ==
nullptr) or (pMEOT->energyLoss() ==
nullptr)) {
806 trajectory.
trackStates().back()->setTrackParameters(
nullptr);
809 std::unique_ptr<const TrackParameters> firstscatpar;
810 std::unique_ptr<const TrackParameters> lastscatpar;
813 double newqoverpid = 0;
816 const double de = std::abs(calomeots[1].energyLoss()->deltaE());
817 const double sigmade = std::abs(calomeots[1].energyLoss()->sigmaDeltaE());
819 const double pbefore = std::abs(1 / firstidpar->parameters()[
Trk::qOverP]);
820 const double pafter = std::abs(1 / tp_closestmuon->parameters()[
Trk::qOverP]);
821 const double elosspull = (pbefore - pafter - de) / sigmade;
823 if (std::abs(elosspull) > 10) {
824 if (elosspull > 10) {
825 newqoverpid = 1 / (de + pafter + 10 * sigmade);
827 newqoverpid = 1 / (de + pafter - 10 * sigmade);
830 if (tp_closestmuon->parameters()[
Trk::qOverP] * newqoverpid < 0) {
834 const AmgVector(5) & newpar = firstidpar->parameters();
835 firstidpar=firstidpar->associatedSurface().createUniqueTrackParameters(
836 newpar[0], newpar[1], newpar[2], newpar[3], newqoverpid, std::nullopt
844 if (lastidpar ==
nullptr) {
850 *(firstismuon ? tp_closestmuon.get() : lastidpar.get()),
851 calomeots[0].associatedSurface(),
858 if (firstscatpar ==
nullptr) {
864 *(firstismuon ? firstidpar : tp_closestmuon),
865 calomeots[2].associatedSurface(),
872 if (lastscatpar ==
nullptr) {
876 std::optional<TransportJacobian> jac1;
877 std::optional<TransportJacobian> jac2;
878 std::unique_ptr<const TrackParameters> elosspar;
880 double firstscatphi = 0;
881 double secondscatphi = 0;
882 double firstscattheta = 0;
883 double secondscattheta = 0;
884 double muonscatphi = 0;
885 double muonscattheta = 0;
887 const TrackParameters *idscatpar = !firstismuon ? firstscatpar.get() : lastscatpar.get();
888 const TrackParameters *muonscatpar = firstismuon ? firstscatpar.get() : lastscatpar.get();
890 newqoverpid = idscatpar->parameters()[
Trk::qOverP];
892 const Amg::Vector3D calosegment = lastscatpar->position() - firstscatpar->position();
895 muonscattheta = calosegment.theta() - muonscatpar->parameters()[
Trk::theta];
898 for (
int i = 0;
i < 2;
i++) {
899 std::unique_ptr<const TrackParameters> tmpelosspar;
900 AmgVector(5) params1 = muonscatpar->parameters();
909 params1[0], params1[1], params1[2], params1[3], params1[4], std::nullopt
930 calomeots[1].associatedSurface(),
936 if ((tmpelosspar ==
nullptr) || (jac1 == std::nullopt)) {
945 const AmgVector(5) & newpars = tmpelosspar->parameters();
946 const std::unique_ptr<const TrackParameters> elosspar2(tmpelosspar->associatedSurface().createUniqueTrackParameters(
947 newpars[0], newpars[1], newpars[2], newpars[3], newqoverpid, std::nullopt
951 elosspar = std::move(tmpelosspar);
954 std::unique_ptr<const TrackParameters> scat2(
m_propagator->propagateParameters(
958 calomeots[0].associatedSurface() :
959 calomeots[2].associatedSurface(),
972 calomeots[0].associatedSurface() :
973 calomeots[2].associatedSurface(),
981 if ((scat2 ==
nullptr) || (jac2 == std::nullopt)) {
986 for (
int j = 0; j < 5; j++) {
987 for (
int k = 0;
k < 5;
k++) {
989 for (
int l = 0;
l < 5;
l++) {
990 jac3[j][
k] += (*jac2) (j,
l) * (*jac1) (
l,
k);
999 jac4(0, 0) = jac3[0][2];
1000 jac4(1, 1) = jac3[1][3];
1001 jac4(0, 1) = jac3[0][3];
1002 jac4(1, 0) = jac3[1][2];
1004 jac4 = jac4.inverse();
1016 discsurf =
static_cast<const Trk::DiscSurface *
>(&scat2->associatedSurface());
1018 if (cylsurf !=
nullptr) {
1022 if (discsurf !=
nullptr) {
1026 double dphi = jac4(0, 0) * dloc1 + jac4(0, 1) * dloc2;
1027 double dtheta = jac4(1, 0) * dloc1 + jac4(1, 1) * dloc2;
1033 muonscatphi += dphi;
1034 muonscattheta += dtheta;
1037 const double idscattheta = idscatpar->parameters()[
Trk::theta] - (scat2->parameters()[
Trk::theta] + dtheta);
1040 firstscatphi = muonscatphi;
1041 secondscatphi = idscatphi;
1042 firstscattheta = muonscattheta;
1043 secondscattheta = idscattheta;
1045 firstscatphi = -idscatphi;
1046 secondscatphi = -muonscatphi;
1047 firstscattheta = -idscattheta;
1048 secondscattheta = -muonscattheta;
1052 AmgVector(5) params2 = scat2->parameters();
1060 firstscatpar=scat2->associatedSurface().createUniqueTrackParameters(
1061 params2[0], params2[1], params2[2], params2[3], params2[4], std::nullopt
1063 idscatpar = firstscatpar.get();
1071 if (startPar !=
nullptr) {
1080 rot.col(2) = trackdir;
1083 trans.linear().matrix() << rot;
1084 trans.translation() << startPar->position() - .1 * trackdir;
1095 if (curvlinpar !=
nullptr) {
1096 startPar.reset(curvlinpar);
1100 firstscatpar = std::move(scat2);
1104 std::unique_ptr<GXFMaterialEffects> firstscatmeff = std::make_unique<GXFMaterialEffects>(calomeots[0]);
1105 std::unique_ptr<GXFMaterialEffects> elossmeff = std::make_unique<GXFMaterialEffects>(calomeots[1]);
1106 std::unique_ptr<GXFMaterialEffects> secondscatmeff = std::make_unique<GXFMaterialEffects>(calomeots[2]);
1108 const double pull1 = std::abs(firstscatphi / firstscatmeff->sigmaDeltaPhi());
1109 const double pull2 = std::abs(secondscatphi / secondscatmeff->sigmaDeltaPhi());
1112 for (
auto &
i : tmp_matvec) {
1117 firstscatmeff->setScatteringAngles(firstscatphi, firstscattheta);
1118 secondscatmeff->setScatteringAngles(secondscatphi, secondscattheta);
1121 elossmeff->setdelta_p(1000 * (lastscatpar->parameters()[
Trk::qOverP] - newqoverpid));
1123 elossmeff->setdelta_p(1000 * (newqoverpid - firstscatpar->parameters()[
Trk::qOverP]));
1126 elossmeff->setSigmaDeltaE(calomeots[1].energyLoss()->sigmaDeltaE());
1128 trajectory.
addMaterialState(std::make_unique<GXFTrackState>(std::move(firstscatmeff), std::move(firstscatpar)), -1);
1129 trajectory.
addMaterialState(std::make_unique<GXFTrackState>(std::move(elossmeff), std::move(elosspar)), -1);
1130 trajectory.
addMaterialState(std::make_unique<GXFTrackState>(std::move(secondscatmeff), std::move(lastscatpar)), -1);
1133 for (
auto &
i : tmp_matvec) {
1140 if (startPar ==
nullptr) {
1145 (pull1 > 5 || pull2 > 5) &&
1151 bool largegap =
false;
1152 double previousz = 0;
1154 for (itStates2 = beginStates2; itStates2 != endState2; ++itStates2) {
1159 if (meff !=
nullptr) {
1165 if ( mefot and mefot->energyLoss() and
1166 std::abs(mefot->energyLoss()->deltaE()) > 250 &&
1167 mefot->energyLoss()->sigmaDeltaE() < 1.e-9
1180 !(itStates2 == beginStates2 || itStates2 == beginStates2 + 1) &&
1189 itStates2 == endState2 - 1 &&
1192 std::abs(tpar->
position().z()) < 13000
1194 std::unique_ptr<const TrackParameters> pseudopar(tpar->
clone());
1198 std::unique_ptr<const PseudoMeasurementOnTrack> newpseudo = std::make_unique<const PseudoMeasurementOnTrack>(
1201 pseudopar->associatedSurface()
1204 std::unique_ptr<GXFTrackState> pseudostate = std::make_unique<GXFTrackState>(std::move(newpseudo), std::move(pseudopar));
1211 pseudostate->setMeasurementErrors(
errors);
1218 std::abs(trajectory.
trackStates().back()->position().z()) > 20000 &&
1219 std::abs(previousz) < 12000
1225 previousz = trajectory.
trackStates().back()->position().z();
1243 const EventContext& ctx,
1245 const Track & intrk1,
1246 const Track & intrk2,
1248 std::vector<MaterialEffectsOnTrack> & calomeots
1250 ATH_MSG_DEBUG(
"--> entering GlobalChi2Fitter::backupCombinationStrategy");
1253 const Track *indettrack = firstismuon ? &intrk2 : &intrk1;
1266 if (pParametersVector->size() > 1)
1267 firstidpar = (*pParametersVector)[1];
1269 firstidpar = pParametersVector->back();
1271 std::unique_ptr<const TrackParameters> lastidpar =
nullptr;
1272 if ((firstidpar !=
nullptr) && (cache.
m_caloEntrance !=
nullptr))
1276 if (lastidpar ==
nullptr) {
1277 lastidpar.reset(pParametersVector->back()->clone());
1280 std::unique_ptr < const TrackParameters > firstscatpar;
1281 std::unique_ptr < const TrackParameters > lastscatpar;
1282 std::unique_ptr < const TrackParameters > elosspar;
1287 lastidpar->position(),
1288 lastidpar->momentum(),
1290 lastidpar->position()
1297 calomeots[0].associatedSurface(),
1303 if (!firstscatpar) {
1307 const std::unique_ptr<const TrackParameters> tmppar(
1311 calomeots[1].associatedSurface(),
1324 const double oldp = std::abs(1 / tmppar->parameters()[
Trk::qOverP]);
1325 const double de = std::abs(calomeots[1].energyLoss()->deltaE());
1327 double newp2 = oldp * oldp - 2 * de * std::sqrt(
mass *
mass + oldp * oldp) + de * de;
1333 const double newqoverp =
sign / std::sqrt(newp2);
1338 tmppar->associatedSurface().createUniqueTrackParameters(
1345 calomeots[2].associatedSurface(),
1358 calomeots[2].associatedSurface(),
1372 calomeots[1].associatedSurface(),
1383 const double sign = (elosspar->parameters()[
Trk::qOverP] < 0) ? -1 : 1;
1384 const double newqoverp =
sign /
1385 (1. / std::abs(elosspar->parameters()[
Trk::qOverP]) +
1386 std::abs(calomeots[1].energyLoss()->deltaE()));
1390 std::unique_ptr<const TrackParameters>
const tmppar(
1391 elosspar->associatedSurface().createUniqueTrackParameters(
1399 calomeots[0].associatedSurface(),
1406 if (!firstscatpar) {
1411 for (; itStates != endState; ++itStates) {
1419 if ((*itStates)->materialEffectsOnTrack() !=
nullptr) {
1432 std::unique_ptr<GXFMaterialEffects> firstscatmeff = std::make_unique<GXFMaterialEffects>(calomeots[0]);
1433 std::unique_ptr<GXFMaterialEffects> elossmeff = std::make_unique<GXFMaterialEffects>(calomeots[1]);
1434 std::unique_ptr<GXFMaterialEffects> secondscatmeff = std::make_unique<GXFMaterialEffects>(calomeots[2]);
1438 sigmadp = calomeots[1].energyLoss()->sigmaDeltaE();
1439 elossmeff->setSigmaDeltaE(sigmadp);
1442 elossmeff->setdelta_p(
dp);
1444 trajectory.
addMaterialState(std::make_unique<GXFTrackState>(std::move(firstscatmeff), std::move(firstscatpar)), -1);
1445 trajectory.
addMaterialState(std::make_unique<GXFTrackState>(std::move(elossmeff), std::move(elosspar)), -1);
1446 trajectory.
addMaterialState(std::make_unique<GXFTrackState>(std::move(secondscatmeff), std::move(lastscatpar)), -1);
1449 const Surface *triggersurf1 =
nullptr;
1450 const Surface *triggersurf2 =
nullptr;
1454 bool seenmdt =
false;
1455 bool mdtbetweenphihits =
false;
1459 itStates2 = (!firstismuon ? beginStates2 : endState - 1);
1460 itStates2 != (!firstismuon ? endState2 : beginStates - 1);
1461 (!firstismuon ? ++itStates2 : --itStates2)
1464 ((*itStates2)->measurementOnTrack() ==
nullptr) ||
1469 const auto *
const pMeasurement = (*itStates2)->measurementOnTrack();
1470 const Surface *surf = &pMeasurement->associatedSurface();
1474 if (isCompetingRIOsOnTrack) {
1476 rot = &crot->rioOnTrack(0);
1479 rot =
static_cast<const RIO_OnTrack *
>(pMeasurement);
1486 (rot !=
nullptr) && (
1493 const double dotprod1 = measdir.dot(
Amg::Vector3D(0, 0, 1));
1496 const bool measphi = std::abs(dotprod1) <= .5 && std::abs(dotprod2) <= .5;
1500 (*itStates2)->trackParameters() !=
nullptr ?
1501 (*itStates2)->trackParameters()->position() :
1503 if (triggersurf1 !=
nullptr) {
1504 triggerpos2 = thispos;
1505 triggersurf2 = surf;
1507 mdtbetweenphihits =
true;
1510 triggerpos1 = thispos;
1511 triggersurf1 = surf;
1517 double mdttrig1 = 999999;
1518 double mdttrig2 = 999999;
1519 const Surface *mdtsurf1 =
nullptr;
1520 const Surface *mdtsurf2 =
nullptr;
1523 itStates2 = (!firstismuon ? beginStates2 : endState - 1);
1524 itStates2 != (!firstismuon ? endState2 : beginStates - 1);
1525 (!firstismuon ? ++itStates2 : --itStates2)
1527 const Surface *surf =
nullptr;
1529 ((*itStates2)->measurementOnTrack() !=
nullptr) &&
1532 surf = &(*itStates2)->measurementOnTrack()->associatedSurface();
1535 if (surf ==
nullptr) {
1538 const auto *
const pThisMeasurement = (*itStates2)->measurementOnTrack();
1542 if (isCompetingRioOnTrack) {
1544 rot = &crot->rioOnTrack(0);
1547 rot =
static_cast<const RIO_OnTrack *
>(pThisMeasurement);
1553 (*itStates2)->trackParameters() !=
nullptr ?
1554 (*itStates2)->trackParameters()->position() :
1555 pThisMeasurement->globalPosition();
1556 if (triggerpos1.mag() > 1 && (globpos - triggerpos1).
mag() < mdttrig1) {
1557 mdttrig1 = (globpos - triggerpos1).
mag();
1560 if (triggerpos2.mag() > 1 && (globpos - triggerpos2).
mag() < mdttrig2) {
1561 mdttrig2 = (globpos - triggerpos2).
mag();
1568 std::vector<GXFTrackState *> outlierstates;
1569 std::vector<GXFTrackState *> outlierstates2;
1571 outlierstates.reserve(10);
1572 outlierstates2.reserve(10);
1574 std::unique_ptr<PseudoMeasurementOnTrack> newpseudo;
1576 for (itStates2 = beginStates2; itStates2 != endState2; ++itStates2) {
1577 const auto *
const pMeasurement{(*itStates2)->measurementOnTrack()};
1579 const bool isStraightLine =
1580 pMeasurement !=
nullptr ?
1587 (newpseudo ==
nullptr) && (
1588 (itStates2 == beginStates2 || itStates2 == beginStates2 + 1) &&
1589 std::abs(pMeasurement->globalPosition().z()) < 10000
1592 std::unique_ptr<const TrackParameters> par2;
1593 if (((*itStates2)->trackParameters() !=
nullptr) && nphi > 99) {
1594 par2.reset((*itStates2)->trackParameters()->clone());
1606 if (par2 ==
nullptr) {
1612 newpseudo = std::make_unique<PseudoMeasurementOnTrack>(
1615 par2->associatedSurface()
1618 std::unique_ptr<GXFTrackState> firstpseudo = std::make_unique<GXFTrackState>(std::move(newpseudo), std::move(par2));
1625 firstpseudo->setMeasurementErrors(
errors);
1626 firstpseudostate = firstpseudo.get();
1632 if (isPseudo && !firstismuon) {
1636 if ((**itStates2).materialEffectsOnTrack() !=
nullptr) {
1646 ((**itStates2).measurementOnTrack() !=
nullptr) &&
1647 &(**itStates2).measurementOnTrack()->associatedSurface() == triggersurf1 &&
1648 (mdtsurf1 !=
nullptr)
1650 std::unique_ptr<Amg::Transform3D> transf = std::make_unique<Amg::Transform3D>(mdtsurf1->
transform());
1652 transf->translation() << triggerpos1;
1657 std::unique_ptr<const PseudoMeasurementOnTrack> newpseudo = std::make_unique<const PseudoMeasurementOnTrack>(
1661 std::unique_ptr<GXFTrackState> pseudostate1 = std::make_unique<GXFTrackState>(std::move(newpseudo),
nullptr);
1668 pseudostate1->setMeasurementErrors(
errors);
1669 outlierstates2.push_back(pseudostate1.get());
1674 ((**itStates2).measurementOnTrack() !=
nullptr) &&
1675 &(**itStates2).measurementOnTrack()->associatedSurface() == triggersurf2 &&
1676 mdtbetweenphihits &&
1677 (mdtsurf2 !=
nullptr)
1679 std::unique_ptr<Amg::Transform3D> transf = std::make_unique<Amg::Transform3D>(mdtsurf2->
transform());
1680 transf->translation() << triggerpos2;
1685 std::unique_ptr<const PseudoMeasurementOnTrack> newpseudo = std::make_unique<const PseudoMeasurementOnTrack>(
1689 std::unique_ptr<GXFTrackState> pseudostate2 = std::make_unique<GXFTrackState>(std::move(newpseudo),
nullptr);
1696 pseudostate2->setMeasurementErrors(
errors);
1698 outlierstates2.push_back(pseudostate2.get());
1714 outlierstates.push_back(trajectory.
trackStates().back().get());
1730 const std::unique_ptr<Trk::Track> tmp_track(
1731 myfit(ctx, cache, trajectory, *startpar2,
false,
muon));
1738 std::abs(trajectory.
residuals().tail<1>()(0) / trajectory.
errors().tail<1>()(0)) > 10
1744 if (firstpseudostate !=
nullptr) {
1749 std::unique_ptr<const PseudoMeasurementOnTrack> newpseudo = std::make_unique<const PseudoMeasurementOnTrack>(
1758 for (
int j = 0; j < (
int) trajectory.
trackStates().size(); j++) {
1759 for (
const auto &
i : outlierstates2) {
1765 for (
const auto &
i : outlierstates) {
1773 itStates = (firstismuon ? beginStates2 : endState - 1);
1774 itStates != (firstismuon ? endState2 : beginStates - 1);
1775 (firstismuon ? ++itStates : --itStates)
1781 makeProtoState(cache, trajectory, *itStates, (firstismuon ? -1 : 0));
1794 std::unique_ptr<Track>
1796 const EventContext& ctx,
1797 const Track& inputTrack,
1801 ATH_MSG_DEBUG(
"--> entering GlobalChi2Fitter::fit(Track,)");
1814 return std::unique_ptr<Track>(
1815 fitIm(ctx, cache, inputTrack, runOutlier, matEffects));
1820 const Track &inputTrack,
1825 const EventContext& ctx = Gaudi::Hive::currentContext();
1829 delete alignCache.m_derivMatrix;
1830 alignCache.m_derivMatrix =
nullptr;
1832 delete alignCache.m_fullCovarianceMatrix;
1833 alignCache.m_fullCovarianceMatrix =
nullptr;
1834 alignCache.m_iterationsOfLastFit = 0;
1837 fitIm(ctx, cache, inputTrack, runOutlier, matEffects);
1838 if(newTrack !=
nullptr){
1843 alignCache.m_iterationsOfLastFit = cache.
m_lastiter;
1850 const EventContext& ctx,
1852 const Track& inputTrack,
1857 ATH_MSG_DEBUG(
"--> entering GlobalChi2Fitter::fit(Track,,)");
1873 ATH_MSG_WARNING(
"Track without track parameters, cannot perform fit");
1881 if (minpar ==
nullptr) {
1899 const Surface *firsthitsurf =
nullptr;
1900 const Surface *lasthitsurf =
nullptr;
1902 bool hasmuon =
false;
1903 bool iscombined =
false;
1904 bool seenphimeas =
false;
1908 for (; itStates != endState; ++itStates) {
1909 const auto *
const pMeasurement = (**itStates).measurementOnTrack();
1911 (pMeasurement ==
nullptr) &&
1912 ((**itStates).materialEffectsOnTrack() !=
nullptr) &&
1918 if (pMeasurement !=
nullptr) {
1919 const Surface *surf = &pMeasurement->associatedSurface();
1928 if (firsthitsurf ==
nullptr) {
1929 firsthitsurf = surf;
1937 if ((**itStates).trackParameters() !=
nullptr) {
1938 lastidpar = (**itStates).trackParameters();
1939 if (firstidpar ==
nullptr) {
1940 firstidpar = lastidpar;
1946 const double dotprod1 = measdir.dot(
Amg::Vector3D(0, 0, 1));
1948 if (std::abs(dotprod1) < .5 && std::abs(dotprod2) < .5) {
1950 if (std::abs(surf->
center().z()) > 13000) {
1953 if (surf->
center().perp() > 9000 && std::abs(surf->
center().z()) < 13000) {
1965 if (iscombined && seenphimeas && (phiem || phibo)) {
1981 if (firstidpar == lastidpar) {
1982 firstidpar = lastidpar =
nullptr;
1992 (firstidpar !=
nullptr)
2003 const bool tmpsirecal = cache.
m_sirecal;
2004 std::unique_ptr<Track> tmptrack =
nullptr;
2015 tmptrack.reset(
myfit(ctx, cache, trajectory, *minpar,
false, matEffects));
2018 if (tmptrack ==
nullptr) {
2027 bool isbrem =
false;
2029 unsigned int n_brem=0;
2031 for (std::unique_ptr<GXFTrackState> & state : trajectory.
trackStates()) {
2034 if (meff !=
nullptr) {
2040 const double p = 1. / std::abs(layerpars->parameters()[
Trk::qOverP] - .0005 * meff->
delta_p());
2042 std::optional<Amg::Vector2D> locpos(state->associatedSurface().globalToLocal(layerpars->
position()));
2043 const Amg::Vector3D layerNormal(state->associatedSurface().normal(*locpos));
2044 double costracksurf = 1.;
2046 costracksurf = std::abs(layerNormal.dot(layerpars->
momentum().unit()));
2048 const double oldde = meff->
deltaE();
2050 std::unique_ptr<EnergyLoss> eloss;
2051 double sigmascat = 0;
2053 if (matprop !=
nullptr) {
2054 eloss = std::make_unique<EnergyLoss>(
2055 m_elosstool->energyLoss(*matprop,
p, 1. / costracksurf,
2057 sigmascat = std::sqrt(
m_scattool->sigmaSquare(*matprop,
p, 1. / costracksurf, matEffects));
2059 if (eloss !=
nullptr) {
2064 sigmascat = std::sqrt(
m_scattool->sigmaSquare(tmpprop,
p, 1. / costracksurf, matEffects));
2082 }
else if (eloss !=
nullptr) {
2093 refpars[0], refpars[1], refpars[2], refpars[3], refpars[4], std::nullopt
2101 trajectory.
brems().clear();
2103 trajectory.
brems().resize(1);
2104 trajectory.
brems()[0] = bremdp;
2115 std::unique_ptr<Track>
track(
myfit(ctx, cache, trajectory, *minpar, runOutlier, matEffects));
2117 bool pseudoupdated =
false;
2119 if ((
track !=
nullptr) && hasid && hasmuon) {
2120 for (std::unique_ptr<GXFTrackState> & pseudostate : trajectory.
trackStates()) {
2122 (pseudostate ==
nullptr) ||
2124 pseudostate->fitQuality().chiSquared() < 10
2130 const std::unique_ptr<const TrackParameters> updpar(
m_updator->removeFromState(
2132 pseudostate->measurement()->localParameters(),
2133 pseudostate->measurement()->localCovariance()
2136 if (updpar ==
nullptr) {
2143 std::unique_ptr<const PseudoMeasurementOnTrack> newpseudo = std::make_unique<const PseudoMeasurementOnTrack>(
2149 pseudostate->setMeasurement(std::move(newpseudo));
2153 pseudostate->setMeasurementErrors(
errors);
2154 pseudoupdated =
true;
2157 if (pseudoupdated) {
2170 if (
track !=
nullptr) {
2173 track->info().addPatternReco(old_info);
2176 return track.release();
2179 std::unique_ptr<Track>
2181 const EventContext& ctx,
2187 ATH_MSG_DEBUG(
"--> entering GlobalChi2Fitter::fit(PRDS,TP,)");
2190 for (
const auto *prd : prds) {
2191 const Surface & prdsurf = (*prd).detectorElement()->surface((*prd).identify());
2196 plsurf =
static_cast < const PlaneSurface *
>(&prdsurf);
2203 if ((slsurf ==
nullptr) && (plsurf ==
nullptr)) {
2204 ATH_MSG_ERROR(
"Surface is neither PlaneSurface nor StraightLineSurface!");
2209 }
else if (slsurf !=
nullptr) {
2218 }
else if (plsurf !=
nullptr) {
2219 if (param.covariance() !=
nullptr) {
2241 if (rot !=
nullptr) {
2242 rots.push_back(rot);
2246 std::unique_ptr<Track>
track =
2247 fit(ctx, rots, param, runOutlier, matEffects);
2249 for (
const auto *rot : rots) {
2256 std::unique_ptr<Track>
2258 const EventContext& ctx,
2259 const Track& inputTrack,
2264 ATH_MSG_DEBUG(
"--> entering GlobalChi2Fitter::fit(Track,Meas'BaseSet,,)");
2279 if (minpar ==
nullptr) {
2297 for (; itStates != endState; ++itStates) {
2302 (trajectory.
trackStates().back()->materialEffects() !=
nullptr) &&
2303 trajectory.
trackStates().back()->materialEffects()->sigmaDeltaE() > 50.001
2305 trajectory.
trackStates().back()->materialEffects()->setKink(
true);
2311 for (
const auto & measBase : addMeasColl) {
2312 if (measBase ==
nullptr) {
2313 ATH_MSG_WARNING(
"There is an empty MeasurementBase object in the track! Skip this object..");
2321 std::unique_ptr<Track>
track(
myfit(ctx, cache, trajectory, *minpar, runOutlier, matEffects));
2324 if (
track !=
nullptr) {
2325 const double oldqual =
2330 const double newqual =
2331 track->fitQuality()->numberDoF() != 0 ?
2332 track->fitQuality()->chiSquared() /
track->fitQuality()->numberDoF() :
2335 if (
m_extensioncuts && runOutlier && newqual > 2 && newqual > 2 * oldqual) {
2339 track.reset(
nullptr);
2343 if (
track !=
nullptr) {
2352 std::unique_ptr<Track>
2354 const EventContext& ctx,
2360 ATH_MSG_DEBUG(
"--> entering GlobalChi2Fitter::fit(Track,PRDS,)");
2364 for (
const auto *prd : prds) {
2365 const Surface & prdsurf = (*prd).detectorElement()->surface((*prd).identify());
2368 std::unique_ptr<const TrackParameters>
const trackparForCorrect(
2384 rot =
m_ROTcreator->correct(*prd, *trackparForCorrect, ctx);
2387 if (rot !=
nullptr) {
2388 rots.push_back(rot);
2392 std::unique_ptr<Track>
track =
fit(ctx,intrk, rots, runOutlier, matEffects);
2394 for (
const auto *rot : rots) {
2402 const EventContext& ctx,
2408 ATH_MSG_DEBUG(
"--> entering GlobalChi2Fitter::fit(Meas'BaseSet,,)");
2421 for (
const auto *itSet : rots) {
2422 if (itSet ==
nullptr) {
2423 ATH_MSG_WARNING(
"There is an empty MeasurementBase object in the track! Skip this object..");
2429 std::unique_ptr<const TrackParameters> startpar(param.
clone());
2432 matEffects ==
muon &&
2438 myfit(ctx,cache, trajectory, *startpar, runOutlier, matEffects);
2456 std::unique_ptr<const PseudoMeasurementOnTrack> newpseudo = std::make_unique<const PseudoMeasurementOnTrack>(
2462 trajectory.
trackStates().front()->setMeasurement(std::move(newpseudo));
2469 std::unique_ptr<const PseudoMeasurementOnTrack> newpseudo = std::make_unique<const PseudoMeasurementOnTrack>(
2475 trajectory.
trackStates().back()->setMeasurement(std::move(newpseudo));
2482 refpars[0], refpars[1], refpars[2], refpars[3], refpars[4], std::nullopt
2487 myfit(ctx,cache, trajectory, *startpar, runOutlier, matEffects);
2498 refpars[0], refpars[1], refpars[2], refpars[3], refpars[4], std::nullopt
2505 firstpar = trajectory.
trackStates().front()->trackParameters();
2510 std::unique_ptr<const PseudoMeasurementOnTrack> newpseudo = std::make_unique<const PseudoMeasurementOnTrack>(
2516 trajectory.
trackStates().front()->setMeasurement(std::move(newpseudo));
2524 lastpar = trajectory.
trackStates().back()->trackParameters();
2529 std::unique_ptr<const PseudoMeasurementOnTrack> newpseudo = std::make_unique<const PseudoMeasurementOnTrack>(
2535 trajectory.
trackStates().back()->setMeasurement(std::move(newpseudo));
2543 std::unique_ptr<Track>
track;
2545 if (startpar !=
nullptr) {
2546 track.reset(
myfit(ctx,cache, trajectory, *startpar, runOutlier, matEffects));
2549 if (
track !=
nullptr) {
2579 std::unique_ptr<GXFMaterialEffects> newmeff;
2587 newmeff = std::make_unique<GXFMaterialEffects>(*meff);
2591 const double sigmascat = std::sqrt(
m_scattool->sigmaSquare(
2612 newmeff = std::make_unique<GXFMaterialEffects>(newmeot);
2628 ((**trajectory.
trackStates().rbegin()).trackParameters() !=
nullptr)
2630 const double delta_p = 1000 * (
2632 (**trajectory.
trackStates().rbegin()).trackParameters()->
2636 newmeff->setdelta_p(delta_p);
2672 seg =
static_cast<const Segment *
>(measbase);
2682 for (
int i = 0;
i <
imax;
i++) {
2685 std::unique_ptr<GXFTrackState> ptsos = std::make_unique<GXFTrackState>(
2686 std::unique_ptr<const MeasurementBase>(measbase2->
clone()),
2687 std::unique_ptr<const TrackParameters>(newtrackpar !=
nullptr ? newtrackpar->
clone() :
nullptr)
2690 double sinstereo = 0;
2703 bool measphi =
false;
2706 bool rotated =
false;
2743 const auto [covEigenValueSmall, covStereoAngle] = principalComponentAnalysis2x2(
covmat);
2744 errors[0] = std::sqrt(covEigenValueSmall);
2745 sinstereo =
std::sin(covStereoAngle);
2760 const double dotprod1 = measdir.dot(
Amg::Vector3D(0, 0, 1));
2762 if (std::abs(dotprod1) < .5 && std::abs(dotprod2) < .5) {
2770 int param_index = 0;
2777 errors[1] = std::sqrt(
covmat(param_index, param_index));
2782 errors[2] = std::sqrt(
covmat(param_index, param_index));
2787 errors[3] = std::sqrt(
covmat(param_index, param_index));
2792 errors[4] = std::sqrt(
covmat(param_index, param_index));
2813 ptsos->setMeasurementErrors(
errors);
2814 ptsos->setSinStereo(sinstereo);
2815 ptsos->setMeasurementType(hittype);
2816 ptsos->setMeasuresPhi(measphi);
2837 if (tvol ==
nullptr) {
2844 if (confinedLayers !=
nullptr) {
2848 if (
layer !=
nullptr) {
2853 if ((layIndex.
value() == 0) || (
layer->layerMaterialProperties() ==
nullptr)) {
2867 if (disclay !=
nullptr) {
2868 if (disclay->
center().z() < 0) {
2873 }
else if (cyllay !=
nullptr) {
2884 for (
size_t ib = 0 ;
ib < bsurf.size(); ++
ib) {
2885 const Layer *
layer = bsurf[
ib]->surfaceRepresentation().materialLayer();
2887 if (
layer ==
nullptr)
continue;
2891 if ((layIndex.
value() == 0) || (
layer->layerMaterialProperties() ==
nullptr)) {
2903 discsurf =
static_cast<const DiscSurface *
>(&
layer->surfaceRepresentation());
2905 if (discsurf !=
nullptr) {
2907 discsurf->
center().z() < 0 &&
2912 discsurf->
center().z() > 0 &&
2918 (cylsurf !=
nullptr) &&
2924 if ((cylsurf ==
nullptr) && (discsurf ==
nullptr)) {
2931 if (confinedVolumes !=
nullptr) {
2932 for (
const auto & volume : confinedVolumes->
arrayObjects()) {
2933 if (volume !=
nullptr) {
2949 std::optional<std::pair<Amg::Vector3D, double>>
2962 const double *
pos = parforextrap.
position().data();
2972 const double r = (std::abs(currentqoverp) > 1
e-10) ? -sintheta / (currentqoverp * 300. *
field[2]) : 1e6;
2973 const double xc = parforextrap.
position().x() -
r * sinphi;
2974 const double yc = parforextrap.
position().y() +
r * cosphi;
2975 const double phi0 = std::atan2(parforextrap.
position().y() - yc, parforextrap.
position().x() - xc);
2976 const double z0 = parforextrap.
position().z();
2977 const double delta_s = (surf.
center().z() -
z0) / costheta;
2978 const double delta_phi = delta_s * sintheta /
r;
2989 const double costracksurf = std::abs(costheta);
2991 return std::make_pair(
intersect, costracksurf);
2994 std::optional<std::pair<Amg::Vector3D, double>>
3009 const double *
pos = parforextrap.
position().data();
3017 const double r = (std::abs(currentqoverp) > 1
e-10) ? -sintheta / (currentqoverp * 300. *
field[2]) : 1e6;
3018 const double xc = parforextrap.
position().x() -
r * sinphi;
3019 const double yc = parforextrap.
position().y() +
r * cosphi;
3020 const double phi0 = std::atan2(parforextrap.
position().y() - yc, parforextrap.
position().x() - xc);
3021 const double z0 = parforextrap.
position().z();
3022 const double d = xc * xc + yc * yc;
3023 const double rcyl = surf.
bounds().
r();
3024 double mysqrt = ((
r + rcyl) * (
r + rcyl) -
d) * (
d - (
r - rcyl) * (
r - rcyl));
3030 mysqrt = std::sqrt(mysqrt);
3031 double firstterm = xc / 2 + (xc * (rcyl * rcyl -
r *
r)) / (2 *
d);
3032 double secondterm = (mysqrt * yc) / (2 *
d);
3033 const double x1 = firstterm + secondterm;
3034 const double x2 = firstterm - secondterm;
3035 firstterm = yc / 2 + (yc * (rcyl * rcyl -
r *
r)) / (2 *
d);
3036 secondterm = (mysqrt * xc) / (2 *
d);
3037 const double y1 = firstterm - secondterm;
3038 const double y2 = firstterm + secondterm;
3041 const double dist1 = (
x -
x1) * (
x -
x1) + (
y -
y1) * (
y -
y1);
3042 const double dist2 = (
x -
x2) * (
x -
x2) + (
y -
y2) * (
y -
y2);
3044 if (dist1 < dist2) {
3052 const double phi1 = std::atan2(
y - yc,
x - xc);
3055 const double delta_z =
r * deltaphi / tantheta;
3056 const double z =
z0 + delta_z;
3069 const double costracksurf = std::abs(normal.unit().dot(trackdir));
3071 return std::make_pair(
intersect, costracksurf);
3078 std::vector<std::pair<const Layer *, const Layer *>> &
layers,
3083 std::vector<std::unique_ptr<GXFTrackState>> &
states = trajectory.
trackStates();
3084 std::vector<std::unique_ptr<GXFTrackState>> oldstates = std::move(
states);
3095 for (
int i = 0;
i <= indexoffset;
i++) {
3105 for (
int i = indexoffset + 1;
i < (
int) oldstates.size();
i++) {
3106 const double rmeas = oldstates[
i]->position().perp();
3107 const double zmeas = oldstates[
i]->position().z();
3115 while (layerindex < (
int)
layers.size()) {
3117 double costracksurf = 0.0;
3138 if (oldstates[
i]->trackParameters() !=
nullptr) {
3139 const double rlayer = cylsurf->
bounds().
r();
3140 if (std::abs(rmeas - rlayer) < std::abs(parforextrap->
position().perp() - rlayer)) {
3141 parforextrap = oldstates[
i]->trackParameters();
3157 if (cylsurf->
bounds().
r() > rmeas)
break;
3167 if (oldstates[
i]->trackParameters() !=
nullptr) {
3168 const double zlayer = discsurf->
center().z();
3169 if (std::abs(zmeas - zlayer) < std::abs(parforextrap->
position().z() - zlayer)) {
3170 parforextrap = oldstates[
i]->trackParameters();
3181 if (std::abs(discsurf->
center().z()) > std::abs(zmeas))
break;
3183 throw std::logic_error(
"Unhandled surface.");
3191 if (matprop ==
nullptr) {
3204 const double actualx0 =
X0 / costracksurf;
3205 const double de = -std::abs(
3209 (
m_p != 0.0 ? std::abs(
m_p) : std::abs(1. / currentqoverp)),
3212 const double sigmascat = std::sqrt(
m_scattool->sigmaSquare(
3214 (
m_p != 0.0 ? std::abs(
m_p) : std::abs(1. / currentqoverp)),
3218 std::unique_ptr<GXFMaterialEffects> meff = std::make_unique<GXFMaterialEffects>();
3219 meff->setDeltaE(de);
3220 meff->setScatteringSigmas(sigmascat / sintheta, sigmascat);
3221 meff->setX0(actualx0);
3222 meff->setSurface(&
layer->surfaceRepresentation());
3223 meff->setMaterialProperties(matprop);
3229 std::unique_ptr<EnergyLoss> eloss;
3232 eloss = std::make_unique<EnergyLoss>(
m_elosstool->energyLoss(
3234 (
m_p != 0.0 ? std::abs(
m_p) : std::abs(1. / currentqoverp)),
3239 if (eloss !=
nullptr) {
3240 meff->setSigmaDeltaE(eloss->sigmaDeltaE());
3245 meff->setDeltaE(-5);
3247 meff->setScatteringSigmas(0, 0);
3250 meff->setSigmaDeltaE(50);
3251 if (eloss !=
nullptr) {
3252 meff->setSigmaDeltaEPos(eloss->sigmaPlusDeltaE());
3253 meff->setSigmaDeltaENeg(eloss->sigmaMinusDeltaE());
3258 "X0: " << meff->x0() <<
" qoverp: " << currentqoverp <<
3259 " sigmascat " << meff->sigmaDeltaTheta() <<
" eloss: " << meff->deltaE() <<
3260 " sigma eloss: " << meff->sigmaDeltaE()
3267 std::unique_ptr<GXFTrackState> matstate = std::make_unique<GXFTrackState>(
3269 std::unique_ptr<const TrackParameters>()
3287 std::vector<std::pair<const Layer *, const Layer *>> &
layers,
3288 std::vector<std::pair<const Layer *, const Layer *>> & upstreamlayers,
3289 const std::vector<std::unique_ptr<GXFTrackState>> & oldstates,
3298 upstreamlayers.reserve(5);
3317 const double lastz = laststate->
position().z();
3318 const double lastr = laststate->
position().perp();
3325 const double slope = (tantheta != 0) ? 1 / tantheta : 0;
3331 std::vector < const Layer *>::const_iterator
it;
3332 std::vector < const Layer *>::const_iterator itend;
3350 for (;
it != itend; ++
it) {
3355 if (std::abs((*it)->surfaceRepresentation().center().z()) > std::abs(lastz)) {
3363 const DiscBounds *discbounds =
static_cast<const DiscBounds *
> (&(*it)->surfaceRepresentation().bounds());
3368 if (discbounds->
rMax() < firstr || discbounds->
rMin() > lastr) {
3372 const double rintersect = firstr + ((*it)->surfaceRepresentation().center().z() - firstz) / slope;
3375 rintersect < discbounds->rMin() - 50 ||
3376 rintersect > discbounds->
rMax() + 50
3386 if ((*
it) == endlayer) {
3398 std::abs((*it)->surfaceRepresentation().center().z()) < std::abs(firstz) ||
3401 upstreamlayers.emplace_back((
Layer *)
nullptr, (*it));
3408 (*
it) != startlayer &&
3409 (std::abs((*it)->surfaceRepresentation().center().z()) > std::abs(firstz2) ||
3410 (*it) == startlayer2)
3424 if (barrelcylinder->surfaceRepresentation().bounds().r() > lastr) {
3431 const double zintersect = firstz + (barrelcylinder->surfaceRepresentation().bounds().r() - firstr) * slope;
3433 if (std::abs(zintersect - barrelcylinder->surfaceRepresentation().center().z()) >
3434 ((
const CylinderSurface*)(&barrelcylinder->surfaceRepresentation()))->bounds().halflengthZ() + 50) {
3438 if (barrelcylinder == endlayer) {
3445 if (barrelcylinder->surfaceRepresentation().bounds().r() < firstr ||
3446 barrelcylinder == startlayer) {
3447 upstreamlayers.emplace_back(barrelcylinder, (
Layer*)
nullptr);
3450 if (barrelcylinder != startlayer &&
3451 (barrelcylinder->surfaceRepresentation().bounds().r() > firstr2 ||
3452 barrelcylinder == startlayer2)) {
3453 layers.emplace_back(barrelcylinder, (
Layer*)
nullptr);
3462 std::sort(upstreamlayers.begin(), upstreamlayers.end(),
GXF::LayerSort());
3466 const EventContext& ctx,
3477 if (!caloEntranceIsValid) {
3504 addMaterial(ctx, cache, trajectory, refpar2, matEffects);
3520 bool hasmat =
false;
3521 int indexoffset = 0;
3522 int lastmatindex = 0;
3523 std::vector<std::unique_ptr<GXFTrackState>> & oldstates = trajectory.
trackStates();
3537 for (
int i = 0;
i < (
int) oldstates.size();
i++) {
3538 if (oldstates[
i]->materialEffects() !=
nullptr) {
3547 if (firstsistate ==
nullptr) {
3548 if (oldstates[
i]->trackParameters() ==
nullptr) {
3549 std::unique_ptr<const TrackParameters> tmppar(
m_propagator->propagateParameters(
3552 oldstates[
i]->associatedSurface(),
3559 if (tmppar ==
nullptr)
return;
3561 oldstates[
i]->setTrackParameters(std::move(tmppar));
3563 firstsistate = oldstates[
i].get();
3565 lastsistate = oldstates[
i].get();
3573 if (lastsistate ==
nullptr) {
3574 throw std::logic_error(
"No track state");
3584 std::unique_ptr<const TrackParameters> tmppar(
m_propagator->propagateParameters(
3593 if (tmppar ==
nullptr)
return;
3605 indexoffset = lastmatindex;
3620 std::vector<std::pair<const Layer *, const Layer *>>
layers;
3621 std::vector<std::pair<const Layer *, const Layer *>> & upstreamlayers = trajectory.
upstreamMaterialLayers();
3635 const EventContext& ctx,
3641 if (refpar2 ==
nullptr) {
3651 std::vector<std::unique_ptr<GXFTrackState>> &
states = trajectory.
trackStates();
3652 std::vector<std::unique_ptr<GXFTrackState>> matstates;
3653 std::unique_ptr< const std::vector < const TrackStateOnSurface *>,
3654 void (*)(
const std::vector<const TrackStateOnSurface *> *) >
3655 matvec(
nullptr,&Trk::GlobalChi2Fitter::Cache::objVectorDeleter<TrackStateOnSurface>);
3656 bool matvec_used=
false;
3657 std::unique_ptr<TrackParameters> startmatpar1;
3658 std::unique_ptr<TrackParameters> startmatpar2;
3669 int npseudomuon1 = 0;
3670 int npseudomuon2 = 0;
3672 for (
auto & state :
states) {
3678 if (firstidhit ==
nullptr) {
3687 if (firsthit ==
nullptr) {
3688 firsthit = state->measurement();
3690 if (
tp ==
nullptr) {
3700 if (
tp ==
nullptr) {
3704 state->setTrackParameters(std::unique_ptr<const TrackParameters>(
tp));
3710 lasthit = state->measurement();
3716 if (firstidhit ==
nullptr) {
3717 firstidhit = state->measurement();
3720 if ((firstidpar ==
nullptr) && (
tp !=
nullptr)) {
3724 lastidhit = state->measurement();
3725 if (
tp !=
nullptr) {
3730 if (firstsiliconpar ==
nullptr) {
3731 firstsiliconpar =
tp;
3733 lastsiliconpar =
tp;
3745 if (firstmuonhit ==
nullptr) {
3746 firstmuonhit = state->measurement();
3747 if (
tp !=
nullptr) {
3751 lastmuonhit = state->measurement();
3752 if (
tp !=
nullptr) {
3758 if (meff->
deltaE() == 0) {
3759 if (firstcalopar ==
nullptr) {
3760 firstcalopar = state->trackParameters();
3762 lastcalopar = state->trackParameters();
3764 if (firstmatpar ==
nullptr) {
3765 firstmatpar = state->trackParameters();
3770 std::unique_ptr<TrackParameters> refpar;
3771 AmgVector(5) newpars = refpar2->parameters();
3778 newpars[0], newpars[1], newpars[2], newpars[3], newpars[4], std::nullopt
3781 if (firstmatpar !=
nullptr) {
3786 if ((startmatpar1 ==
nullptr) || ((firstidhit !=
nullptr) && (firstmuonhit !=
nullptr))) {
3790 const double mass = trajectory.
mass();
3792 const AmgVector(5) & newpars = startmatpar2->parameters();
3793 const double oldp = std::abs(1 / newpars[
Trk::qOverP]);
3796 startmatpar2 = startmatpar2->associatedSurface().createUniqueTrackParameters(
3797 newpars[0], newpars[1], newpars[2], newpars[3],
3798 sign / std::sqrt(oldp * oldp + 2 * 100 *
MeV * std::sqrt(oldp * oldp +
mass *
mass) + 100 *
MeV * 100 *
MeV),
3803 AmgVector(5) newpars = startmatpar1->parameters();
3806 startmatpar1 = startmatpar1->associatedSurface().createUniqueTrackParameters(
3807 newpars[0], newpars[1], newpars[2], newpars[3], newpars[4], std::nullopt
3810 newpars = startmatpar2->parameters();
3813 startmatpar2 = startmatpar2->associatedSurface().createUniqueTrackParameters(
3814 newpars[0], newpars[1], newpars[2], newpars[3], newpars[4], std::nullopt
3822 refpar->momentum().unit()
3828 ATH_MSG_DEBUG(
"Obtaining upstream layers from Extrapolator");
3831 std::unique_ptr<const TrackParameters> tmppar;
3833 if (firstmuonhit !=
nullptr) {
3835 if (caloEntranceIsValid) {
3842 if (tmppar !=
nullptr) {
3843 destsurf = &tmppar->associatedSurface();
3848 if (matvec_used) cache.
m_matTempStore.push_back( std::move(matvec) );
3853 false, matEffects) );
3856 if (matvec && !matvec->empty()) {
3857 for (
int i = (
int)matvec->size() - 1;
i > -1;
i--) {
3862 std::unique_ptr<GXFMaterialEffects> meff = std::make_unique<GXFMaterialEffects>(*meot);
3863 const TrackParameters * newpars = (*matvec)[
i]->trackParameters() !=
nullptr ? (*matvec)[
i]->trackParameters()->
clone() :
nullptr;
3864 meff->setSigmaDeltaE(0);
3865 matstates.push_back(std::make_unique<GXFTrackState>(
3867 std::unique_ptr<const TrackParameters>(newpars)
3880 refpar->momentum().unit()
3886 ATH_MSG_DEBUG(
"Obtaining downstream ID layers from Extrapolator");
3888 std::unique_ptr<const TrackParameters> tmppar;
3889 std::unique_ptr<Surface> calosurf;
3890 if (firstmuonhit !=
nullptr) {
3892 if (caloEntranceIsValid) {
3900 if (tmppar !=
nullptr) {
3904 cylcalosurf =
static_cast<const CylinderSurface *
>(&tmppar->associatedSurface());
3909 disccalosurf =
static_cast<const DiscSurface *
>(&tmppar->associatedSurface());
3911 if (cylcalosurf !=
nullptr) {
3914 const double radius = cylbounds.
r();
3916 calosurf = std::make_unique<CylinderSurface>(trans,
radius - 1, hlength);
3917 }
else if (disccalosurf !=
nullptr) {
3918 const double newz = (
3919 disccalosurf->
center().z() > 0 ?
3920 disccalosurf->
center().z() - 1 :
3921 disccalosurf->
center().z() + 1
3925 disccalosurf->
center().x(),
3926 disccalosurf->
center().y(),
3931 trans.translation() << newpos;
3934 const double rmin = discbounds->
rMin();
3935 const double rmax = discbounds->
rMax();
3936 calosurf = std::make_unique<DiscSurface>(trans, rmin, rmax);
3938 destsurf = calosurf.release();
3942 if (matvec_used) cache.
m_matTempStore.push_back( std::move(matvec) );
3944 ctx, *startmatpar2, *destsurf,
alongMomentum,
false, matEffects));
3945 matvec_used =
false;
3947 if (matvec && !matvec->empty()) {
3948 for (
const auto &
i : *matvec) {
3954 std::unique_ptr<GXFMaterialEffects> meff = std::make_unique<GXFMaterialEffects>(*meot);
3960 meff->setDeltaE(-5);
3963 meff->setScatteringSigmas(0, 0);
3966 meff->setSigmaDeltaE(50);
3969 const TrackParameters * newparams =
i->trackParameters() !=
nullptr ?
i->trackParameters()->clone() :
nullptr;
3971 matstates.push_back(std::make_unique<GXFTrackState>(
3973 std::unique_ptr<const TrackParameters>(newparams)
3985 if (cache.
m_calomat && (firstmuonhit !=
nullptr) && (firstidhit !=
nullptr)) {
3988 std::vector<MaterialEffectsOnTrack> calomeots =
m_calotool->extrapolationSurfacesAndEffects(
3997 if (calomeots.empty()) {
4002 std::unique_ptr<const TrackParameters> prevtrackpars =
unique_clone(lastidpar);
4003 if (lasthit == lastmuonhit) {
4004 for (
int i = 0;
i < (
int) calomeots.size();
i++) {
4007 std::unique_ptr<const TrackParameters> layerpar(
m_propagator->propagateParameters(
4010 calomeots[
i].associatedSurface(),
4017 if (layerpar ==
nullptr) {
4022 std::unique_ptr<GXFMaterialEffects> meff = std::make_unique<GXFMaterialEffects>(calomeots[
i]);
4025 lastcalopar = layerpar.get();
4029 const double qoverp = layerpar->parameters()[
Trk::qOverP];
4030 double qoverpbrem = 0;
4034 (firstmuonpar !=
nullptr) &&
4035 std::abs(firstmuonpar->parameters()[
Trk::qOverP]) > 1.e-9
4037 qoverpbrem = firstmuonpar->parameters()[
Trk::qOverP];
4039 const double sign = (qoverp > 0) ? 1 : -1;
4040 qoverpbrem =
sign / (1 / std::abs(qoverp) - std::abs(calomeots[
i].energyLoss()->deltaE()));
4043 const AmgVector(5) & newpar = layerpar->parameters();
4045 layerpar = layerpar->associatedSurface().createUniqueTrackParameters(
4046 newpar[0], newpar[1], newpar[2], newpar[3], qoverpbrem, std::nullopt
4048 meff->setdelta_p(1000 * (qoverpbrem - qoverp));
4051 matstates.push_back(std::make_unique<GXFTrackState>(
4053 std::unique_ptr<const TrackParameters>(layerpar !=
nullptr ? layerpar->clone() :
nullptr)
4055 prevtrackpars = std::move(layerpar);
4060 firsthit == firstmuonhit &&
4064 for (
int i = 0;
i < (
int) calomeots.size();
i++) {
4066 std::unique_ptr<const TrackParameters> layerpar(
m_propagator->propagateParameters(
4069 calomeots[
i].associatedSurface(),
4076 if (layerpar ==
nullptr) {
4081 std::unique_ptr<GXFMaterialEffects> meff = std::make_unique<GXFMaterialEffects>(calomeots[
i]);
4090 const double qoverpbrem = layerpar->parameters()[
Trk::qOverP];
4095 (lastmuonpar !=
nullptr) &&
4096 std::abs(lastmuonpar->parameters()[
Trk::qOverP]) > 1.e-9
4100 const double sign = (qoverpbrem > 0) ? 1 : -1;
4101 qoverp =
sign / (1 / std::abs(qoverpbrem) + std::abs(calomeots[
i].energyLoss()->deltaE()));
4104 meff->setdelta_p(1000 * (qoverpbrem - qoverp));
4105 const AmgVector(5) & newpar = layerpar->parameters();
4107 prevtrackpars = layerpar->associatedSurface().createUniqueTrackParameters(
4108 newpar[0], newpar[1], newpar[2], newpar[3], qoverp, std::nullopt
4112 matstates.insert(matstates.begin(), std::make_unique<GXFTrackState>(std::move(meff), std::move(layerpar)));
4117 if (lasthit == lastmuonhit && cache.
m_extmat) {
4118 std::unique_ptr<const Trk::TrackParameters> muonpar1;
4120 if (lastcalopar !=
nullptr) {
4122 if (msEntranceIsValid) {
4130 if (muonpar1 !=
nullptr) {
4138 rot.col(2) = trackdir;
4140 trans.linear().matrix() << rot;
4141 trans.translation() << muonpar1->
position() - .1 * trackdir;
4144 std::unique_ptr<const TrackParameters> curvlinpar(
m_extrapolator->extrapolateDirectly(
4152 if (curvlinpar !=
nullptr) {
4153 muonpar1 = std::move(curvlinpar);
4157 muonpar1 = std::unique_ptr<TrackParameters>(lastcalopar->
clone());
4161 muonpar1 = std::unique_ptr<TrackParameters>(refpar->clone());
4166 if (muonpar1 !=
nullptr) {
4176 0) and (firstmuonhit !=
nullptr)) {
4188 std::abs(distsol.
first()) < std::abs(distsol.
second()) ?
4195 if (firstmuonpar !=
nullptr) {
4196 AmgVector(5) newpars = firstmuonpar->parameters();
4203 newpars[0], newpars[1], newpars[2], newpars[3], newpars[4], std::nullopt
4206 std::unique_ptr<const TrackParameters> tmppar(
m_propagator->propagateParameters(
4216 if (tmppar !=
nullptr) {
4217 muonpar1 = std::move(tmppar);
4223 ATH_MSG_DEBUG(
"Obtaining downstream layers from Extrapolator");
4224 if (matvec_used) cache.
m_matTempStore.push_back( std::move(matvec) );
4227 states.back()->associatedSurface(),
4231 matvec_used =
false;
4238 if (matvec && !matvec->empty()) {
4239 for (
int j = 0; j < (
int) matvec->size(); j++) {
4245 std::unique_ptr<GXFMaterialEffects> meff = std::make_unique<GXFMaterialEffects>(*meot);
4250 std::abs(meff->deltaE()) > 25 &&
4251 std::abs((*matvec)[j]->trackParameters()->position().z()) < 13000
4256 const TrackParameters * newparams = (*matvec)[j]->trackParameters() !=
nullptr ? (*matvec)[j]->trackParameters()->
clone() :
nullptr;
4258 matstates.push_back(std::make_unique<GXFTrackState>(
4260 std::unique_ptr<const TrackParameters>(newparams)
4270 if (firsthit == firstmuonhit && cache.
m_extmat && (firstcalopar !=
nullptr)) {
4271 std::unique_ptr<const Trk::TrackParameters> muonpar1;
4274 if (msEntranceIsValid) {
4282 if (muonpar1 !=
nullptr) {
4290 rot.col(2) = trackdir;
4292 trans.linear().matrix() << rot;
4293 trans.translation() << muonpar1->
position() - .1 * trackdir;
4296 std::unique_ptr<const TrackParameters> curvlinpar(
m_extrapolator->extrapolateDirectly(
4304 if (curvlinpar !=
nullptr) {
4305 muonpar1 = std::move(curvlinpar);
4309 muonpar1 = std::unique_ptr<const TrackParameters>(firstcalopar->
clone());
4315 if (muonpar1 !=
nullptr) {
4326 ATH_MSG_DEBUG(
"Collecting upstream muon material from extrapolator");
4327 if (matvec_used) cache.
m_matTempStore.push_back( std::move(matvec) );
4330 states[0]->associatedSurface(),
4334 matvec_used =
false;
4336 if (matvec && !matvec->empty()) {
4337 ATH_MSG_DEBUG(
"Retrieved " << matvec->size() <<
" material states");
4339 for (
int j = 0; j < (
int) matvec->size(); j++) {
4342 if (meb !=
nullptr) {
4347 std::unique_ptr<GXFMaterialEffects> meff = std::make_unique<GXFMaterialEffects>(*meot);
4352 std::abs(meff->deltaE()) > 25 &&
4353 std::abs((*matvec)[j]->trackParameters()->position().z()) < 13000
4359 (*matvec)[j]->trackParameters() !=
nullptr
4360 ? (*matvec)[j]->trackParameters()->
clone()
4363 matstates.insert(matstates.begin(), std::make_unique<GXFTrackState>(
4365 std::unique_ptr<const TrackParameters>(tmpparams)
4378 std::vector<std::unique_ptr<GXFTrackState>> & newstates =
states;
4379 std::vector<std::unique_ptr<GXFTrackState>> oldstates = std::move(newstates);
4382 newstates.reserve(oldstates.size() + matstates.size());
4385 int firstlayerno = -1;
4395 bool addlayer =
true;
4397 while (addlayer && layerno < (
int) matstates.size()) {
4399 const TrackParameters *layerpar = matstates[layerno]->trackParameters();
4401 const DistanceSolution distsol = oldstates[
i]->associatedSurface().straightLineDistanceEstimate(
4414 const double trackimpact = std::abs(-refpar->position().x() * sinphi + refpar->position().y() * cosphi);
4416 if (trackimpact > cylinderradius - 5 *
mm) {
4422 if (
i == (
int) oldstates.size() - 1) {
4431 ATH_MSG_DEBUG(
"Material state with excessive scattering, skipping it");
4442 if (firstlayerno < 0) {
4443 firstlayerno = layerno;
4448 if ((layerpar !=
nullptr) && matEffects !=
pion && matEffects !=
muon) {
4450 const Layer *lay =
nullptr;
4452 if (tvol !=
nullptr) {
4468 if (matvec_used) cache.
m_matTempStore.push_back( std::move(matvec) );
4493 std::unique_ptr<const TrackParameters> per(
m_extrapolator->extrapolate(
4494 Gaudi::Hive::currentContext(),param, tmppersf,
oppositeMomentum,
false, matEffects));
4496 if (per ==
nullptr) {
4497 ATH_MSG_DEBUG(
"Cannot make Perigee with starting parameters");
4501 if(std::abs(per->position().z())>5000.) {
4502 ATH_MSG_WARNING(
"Pathological perigee well outside of tracking detector!! Returning nullptr");
4510 const EventContext& ctx,
4517 ATH_MSG_DEBUG(
"--> entering GlobalChi2Fitter::myfit_helper");
4545 if (trajectory.
nDOF() < 0) {
4555 ATH_MSG_WARNING(
"Attempt to apply material corrections with q/p=0, reject track");
4569 std::unique_ptr<const TrackParameters> per =
makePerigee(cache, param, matEffects);
4585 ATH_MSG_WARNING(
"Tracking Geometry Service and/or Material Updator Tool not configured");
4596 addMaterial(ctx, cache, trajectory, per !=
nullptr ? per.get() : ¶m, matEffects);
4599 ctx, cache, trajectory, per !=
nullptr ? per.get() : ¶m, matEffects);
4611 for (
const auto & state : trajectory.
trackStates()) {
4615 state->materialEffects()->deltaE() == 0
4617 scatstate2 = state.get();
4622 scatstate = state.get();
4629 if ((scatstate ==
nullptr) || (scatstate2 ==
nullptr)) {
4630 throw std::logic_error(
"Invalid scatterer");
4637 trajectory.
trackStates()[nstates / 2 - 1]->position() +
4643 std::unique_ptr<const TrackParameters> nearestpar;
4644 double mindist = 99999;
4645 std::vector < GXFTrackState * >mymatvec;
4648 if ((*it).trackParameters() ==
nullptr) {
4654 (*it).trackParameters()->position(),
4655 (*it).trackParameters()->momentum().unit())
4658 const bool insideid = (
4664 (((*it).measurement() !=
nullptr) && insideid) || (
4665 ((*it).materialEffects() !=
nullptr) &&
4667 (*it).materialEffects()->deltaE() == 0 ||
4668 ((*it).materialEffects()->sigmaDeltaPhi() == 0 &&
4670 (*it).materialEffects()->deltaPhi() != 0
4674 const double dist = ((*it).trackParameters()->position() -
vertex).
perp();
4675 if (dist < mindist) {
4683 if (((*it).materialEffects() !=
nullptr) &&
distance > 0) {
4684 mymatvec.push_back(
it.get());
4688 if (nearestpar ==
nullptr) {
4692 for (
auto & state : mymatvec) {
4694 const Surface &matsurf = state->associatedSurface();
4696 nearestpar->position(), nearestpar->momentum().unit());
4704 std::unique_ptr<const TrackParameters> tmppar(
m_propagator->propagateParameters(
4714 if (tmppar ==
nullptr) {
4726 if (tmppar ==
nullptr) {
4736 AmgVector(5) newpars = tmppar->parameters();
4738 if (state->materialEffects()->sigmaDeltaE() > 0) {
4739 newpars[
Trk::qOverP] += .001 * state->materialEffects()->delta_p();
4742 const double de = std::abs(state->materialEffects()->deltaE());
4743 const double oldp = std::abs(1 / newpars[
Trk::qOverP]);
4744 const double newp2 = oldp * oldp - 2 * de * std::sqrt(
mass *
mass + oldp * oldp) + de * de;
4750 nearestpar = tmppar->associatedSurface().createUniqueTrackParameters(
4751 newpars[0], newpars[1], newpars[2], newpars[3], newpars[4], std::nullopt
4755 std::unique_ptr<Trk::TrackParameters> tmpPars(
m_propagator->propagateParameters(
4766 if (tmpPars !=
nullptr) {
4767 per.reset(
static_cast < const Perigee *
>(tmpPars.release()));
4772 const double oldp = 1. / std::abs(per->parameters()[
Trk::qOverP]);
4774 const double newp = std::sqrt(oldp * oldp + 2 * toteloss * std::sqrt(oldp * oldp +
mass *
mass) + toteloss * toteloss);
4778 per = per->associatedSurface().createUniqueTrackParameters(
4783 if (per ==
nullptr) {
4800 }
else if (per ==
nullptr) {
4818 per = per->associatedSurface().createUniqueTrackParameters(
4825 if (per !=
nullptr) {
4835 Eigen::MatrixXd a_inv;
4836 a.resize(nfitpar, nfitpar);
4841 derivPool.setZero();
4843 for (std::unique_ptr<GXFTrackState> & state : trajectory.
trackStates()) {
4844 if (state->materialEffects() !=
nullptr) {
4847 state->setDerivatives(derivPool);
4850 bool doderiv =
true;
4882 const double redchi2 = (trajectory.
nDOF() > 0) ? trajectory.
chi2() / trajectory.
nDOF() : 0;
4883 const double prevredchi2 = (trajectory.
nDOF() > 0) ? trajectory.
prevchi2() / trajectory.
nDOF() : 0;
4892 (redchi2 < prevredchi2 &&
4893 (redchi2 > prevredchi2 - 1 || redchi2 < 2)) ||
4894 nsihits + ntrthits == nhits
4899 if (
it != 1 || nsihits != 0 || trajectory.
nDOF() <= 0 || trajectory.
chi2() / trajectory.
nDOF() <= 3) {
4914 if (ntrtprechits+ntrttubehits) {
4915 phf =
float(ntrtprechits)/
float(ntrtprechits+ntrttubehits);
4917 if (phf<m_minphfcut && it>=3) {
4918 if ((ntrtprechits+ntrttubehits)>=15) {
4923 <<
" | nTRTPrecHits = " << ntrtprechits
4924 <<
" | nTRTTubeHits = " << ntrttubehits
4925 <<
" | nOutliers = "
4945 if (trajectory.
prefit() == 0) {
4949 Eigen::LLT < Eigen::MatrixXd >
const lltOfW(
a);
4950 if (lltOfW.info() == Eigen::Success) {
4954 const int ncols =
a.cols();
4955 Amg::MatrixX const weightInvAMG = Amg::MatrixX::Identity(ncols, ncols);
4956 a_inv = lltOfW.solve(weightInvAMG);
4975 if (traj != &trajectory) {
4980 finaltrajectory = traj;
4988 if (nperpars == 5) {
4989 for (
int i = 0;
i <
a.cols();
i++) {
4990 a_inv(4,
i) *= .001;
4991 a_inv(
i, 4) *= .001;
4995 int scatterPos = nperpars + 2 * nscat;
4996 for (
int bremno = 0; bremno < nbrem; bremno++, scatterPos++) {
4997 for (
int i = 0;
i <
a.cols();
i++) {
4998 a_inv(scatterPos,
i) *= .001;
4999 a_inv(
i, scatterPos) *= .001;
5006 for (
int i = 0;
i < nperparams;
i++) {
5007 for (
int j = 0; j < nperparams; j++) {
5008 (errmat) (j,
i) = a_inv(j,
i);
5013 (errmat) (4, 4) = 1
e-20;
5017 std::unique_ptr<const TrackParameters> measper(
5019 perpars[0], perpars[1], perpars[2], perpars[3], perpars[4], std::move(errmat)
5029 std::unique_ptr<Track>
track =
nullptr;
5031 if (finaltrajectory->
prefit() > 0) {
5032 if (finaltrajectory != &trajectory) {
5034 delete finaltrajectory;
5053 (
track !=
nullptr) && (
5054 track->fitQuality()->numberDoF() != 0 &&
5055 track->fitQuality()->chiSquared() /
track->fitQuality()->numberDoF() >
cut
5058 track.reset(
nullptr);
5062 if (
track ==
nullptr) {
5066 if (finaltrajectory != &trajectory) {
5067 delete finaltrajectory;
5070 return track.release();
5074 const EventContext& ctx,
5075 const Cache & cache,
5079 int & bremno_maxbrempull,
5084 std::vector<std::unique_ptr<GXFTrackState>> &
states = trajectory.
trackStates();
5104 const int nmeas = (
int)
res.size();
5114 const int nDOF = trajectory.
nDOF();
5115 const bool doNewPseudoMeasurements = (
5119 std::abs((trajectory.
prevchi2() - trajectory.
chi2()) / nDOF) < 15 &&
5121 (nperpars == 0 || nidhits > 0)
5132 double maxbrempull = -0.2;
5141 for (
int hitno = 0; hitno < (
int)
states.size(); hitno++) {
5142 std::unique_ptr<GXFTrackState> & state =
states[hitno];
5156 doNewPseudoMeasurements &&
5158 !state->associatedSurface().isFree() &&
5159 !state->isRecalibrated()
5164 std::unique_ptr<const PseudoMeasurementOnTrack> newpseudo = std::make_unique<const PseudoMeasurementOnTrack>(
5170 state->setMeasurement(std::move(newpseudo));
5171 measbase = state->measurement();
5178 double *
errors = state->measurementErrors();
5180 for (
int i = 0;
i < 5;
i++) {
5200 res[measno] = residuals[
i];
5215 double *
errors = state->measurementErrors();
5216 for (
int i = 0;
i < 5;
i++) {
5229 ((trajectory.
prefit() == 0) || state->materialEffects()->deltaE() == 0)
5231 const double deltaPhi = state->materialEffects()->deltaPhi();
5232 const double measDeltaPhi = state->materialEffects()->measuredDeltaPhi();
5233 const double sigma2deltaPhi =
std::pow(state->materialEffects()->sigmaDeltaPhi(), 2);
5234 const double deltaTheta = state->materialEffects()->deltaTheta();
5235 const double sigma2deltaTheta =
std::pow(state->materialEffects()->sigmaDeltaTheta(), 2);
5237 if (trajectory.
prefit() != 1) {
5238 b[nperpars + 2 * scatno] -= (
deltaPhi - measDeltaPhi) / sigma2deltaPhi;
5239 b[nperpars + 2 * scatno + 1] -= deltaTheta / sigma2deltaTheta;
5241 b[nperpars + scatno] -= deltaTheta / sigma2deltaTheta;
5246 deltaTheta * deltaTheta / sigma2deltaTheta
5255 if ((state->materialEffects() !=
nullptr) && state->materialEffects()->sigmaDeltaE() > 0) {
5256 double averagenergyloss = std::abs(state->materialEffects()->deltaE());
5258 const double qoverp = limitInversePValue(qoverpbrem - state->materialEffects()->delta_p());
5259 const double pbrem = 1. / std::abs(qoverpbrem);
5260 const double p = 1. / std::abs(qoverp);
5261 const double mass = .001 * trajectory.
mass();
5263 const double bremEnergy = std::sqrt(pbrem * pbrem +
mass *
mass);
5265 const double resMaterial = .001 * averagenergyloss -
energy + bremEnergy;
5266 res[nmeas - nbrem + bremno] = resMaterial;
5268 const double sigde = state->materialEffects()->sigmaDeltaE();
5269 const double sigdepos = state->materialEffects()->sigmaDeltaEPos();
5270 const double sigdeneg = state->materialEffects()->sigmaDeltaENeg();
5272 double errorMaterial = .001 * state->materialEffects()->sigmaDeltaE();
5273 error[nmeas - nbrem + bremno] = errorMaterial;
5284 if (state->materialEffects()->isKink()) {
5285 maxbrempull = -999999999;
5286 state_maxbrempull =
nullptr;
5292 trajectory.
prefit() == 0 &&
5294 sigde != sigdepos &&
5297 const double elosspull = resMaterial / errorMaterial;
5299 if (trajectory.
mass() > 100) {
5304 if (std::abs(elosspull) > 1) {
5305 if (elosspull < -1) {
5306 state->materialEffects()->setSigmaDeltaE(sigdepos);
5308 state->materialEffects()->setSigmaDeltaE(sigdeneg);
5311 errorMaterial = .001 * state->materialEffects()->sigmaDeltaE();
5312 error[nmeas - nbrem + bremno] = errorMaterial;
5324 !state->materialEffects()->isKink() && (
5325 (
m_fixbrem == -1 && elosspull < maxbrempull) ||
5329 bremno_maxbrempull = bremno;
5330 state_maxbrempull = state.get();
5331 maxbrempull = elosspull;
5340 trajectory.
prefit() == 0 &&
5341 state->materialEffects()->sigmaDeltaPhi() == 0 &&
5342 state->materialEffects()->isMeasuredEloss() &&
5343 resMaterial / (.001 * state->materialEffects()->sigmaDeltaEAve()) > 2.5
5348 std::vector<MaterialEffectsOnTrack> calomeots =
5361 if (calomeots.size() == 3) {
5362 averagenergyloss = std::abs(calomeots[1].energyLoss()->deltaE());
5363 const double newres = .001 * averagenergyloss -
energy + bremEnergy;
5364 const double newerr = .001 * calomeots[1].energyLoss()->sigmaDeltaE();
5366 const double oldPull = resMaterial / errorMaterial;
5367 const double newPull = newres / newerr;
5369 if (std::abs(newPull) < std::abs(oldPull)) {
5370 ATH_MSG_DEBUG(
"Changing from measured to parametrized energy loss");
5372 state->materialEffects()->setEloss(std::unique_ptr<EnergyLoss>(calomeots[1].energyLoss()->
clone()));
5373 state->materialEffects()->setSigmaDeltaE(calomeots[1].energyLoss()->sigmaDeltaE());
5374 res[nmeas - nbrem + bremno] = newres;
5375 error[nmeas - nbrem + bremno] = newerr;
5379 state->materialEffects()->setMeasuredEloss(
false);
5389 for (
int imeas = 0; imeas < nmeas; imeas++) {
5390 if (
error[imeas] == 0) {
5405 const Cache & cache,
5411 const double oldChi2 = trajectory.
prevchi2();
5412 const double newChi2 = trajectory.
chi2();
5417 const double nDOF = trajectory.
nDOF();
5418 const double oldRedChi2 = (nDOF > 0) ? oldChi2 / nDOF : 0;
5419 const double newRedChi2 = (nDOF > 0) ? newChi2 / nDOF : 0;
5422 trajectory.
prefit() > 0 && (
5423 (newRedChi2 < 2 &&
it != 0) ||
5424 (newRedChi2 < oldRedChi2 + .1 && std::abs(newRedChi2 - oldRedChi2) < 1 &&
it != 1)
5437 int miniter = (nsihits != 0 && nsihits + ntrthits == nhits) ? 1 : 2;
5440 if (
it >= miniter && std::abs(oldChi2 - newChi2) < 1) {
5447 const int bremno_maxbrempull,
5453 if (state_maxbrempull ==
nullptr) {
5465 const int nmeas = (
int)
res.size();
5468 const double oldError =
error[nmeas - nbrem + bremno_maxbrempull];
5470 error[nmeas - nbrem + bremno_maxbrempull] = newError;
5473 if (
a.cols() != nFitPars) {
5477 const double errorRatio = oldError / newError;
5478 const double errorReductionRatio = 1 -
std::pow(errorRatio, 2);
5481 for (
int i = 0;
i < nFitPars;
i++) {
5482 if (weightderiv(nmeas - nbrem + bremno_maxbrempull,
i) == 0) {
5486 for (
int j =
i; j < nFitPars; j++) {
5487 const double newaij =
a(
i, j) - errorReductionRatio *
5488 weightderiv(nmeas - nbrem + bremno_maxbrempull,
i) *
5489 weightderiv(nmeas - nbrem + bremno_maxbrempull, j);
5491 a.fillSymmetric(
i, j, newaij);
5493 weightderiv(nmeas - nbrem + bremno_maxbrempull,
i) *= errorRatio;
5502 std::vector<std::unique_ptr<GXFTrackState>> &
states = trajectory.
trackStates();
5515 const int nmeas = (
int) weightderiv.rows();
5517 for (std::unique_ptr<GXFTrackState> & state :
states) {
5521 const auto [scatmin, scatmax] = std::minmax(scatno, nscatupstream);
5522 const auto [bremmin, bremmax] = std::minmax(bremno, nbremupstream);
5529 const double sinStereo =
5531 state->sinStereo() :
5533 const double cosStereo =
5535 std::sqrt(1 -
std::pow(sinStereo, 2)) :
5543 auto getThisDeriv = [sinStereo, cosStereo, &derivatives](
int i,
int j) ->
double {
5544 if (
i == 0 && sinStereo != 0) {
5545 return derivatives(0, j) * cosStereo + sinStereo * derivatives(1, j);
5547 return derivatives(
i, j);
5551 for (
int i = 0;
i < 5;
i++) {
5566 if (
i == 0 && sinStereo != 0) {
5567 weightderiv.row(measno).head(
cols) =
5568 (derivatives.row(0).head(
cols) * cosStereo +
5569 sinStereo * derivatives.row(1).head(
cols)) /
5572 weightderiv.row(measno).head(
cols) = derivatives.row(
i).head(
cols) /
error[measno];
5576 for (
int j = scatmin; j < scatmax; j++) {
5577 if (trajectory.
prefit() == 1) {
5578 const int index = nperparams + j;
5581 const int index = nperparams + 2 * j;
5583 weightderiv(measno,
index + 1) = getThisDeriv(
i,
index + 1) /
error[measno];
5587 for (
int j = bremmin; j < bremmax; j++) {
5588 const int index = j + nperparams + 2 * nscat;
5595 double *
errors = state->measurementErrors();
5596 for (
int i = 0;
i < 5;
i++) {
5603 ((trajectory.
prefit() == 0) || state->materialEffects()->deltaE() == 0)
5608 if ((state->materialEffects() !=
nullptr) && state->materialEffects()->sigmaDeltaE() > 0) {
5610 const double qoverpbrem = limitInversePValue(1000 * state->trackParameters()->parameters()[
Trk::qOverP]);
5611 const double qoverp = limitInversePValue(qoverpbrem - state->materialEffects()->delta_p());
5613 const double mass = .001 * trajectory.
mass();
5615 const auto thisMeasurementIdx{nmeas - nbrem + bremno};
5619 auto multiplier = [] (
double mass,
double qOverP){
5622 const auto qoverpTerm {multiplier(
mass, qoverp) /
error[thisMeasurementIdx]};
5623 const auto qoverpBremTerm {multiplier(
mass, qoverpbrem) /
error[thisMeasurementIdx]};
5626 weightderiv(thisMeasurementIdx, 4) = qoverpBremTerm - qoverpTerm;
5629 const auto bremNoBase = nperparams + 2 * nscat;
5630 if (bremno < nbremupstream) {
5631 weightderiv(thisMeasurementIdx, bremNoBase + bremno) = qoverpTerm;
5632 for (
int bremno2 = bremno + 1; bremno2 < nbremupstream; bremno2++) {
5633 weightderiv(thisMeasurementIdx, bremNoBase + bremno2) = qoverpTerm - qoverpBremTerm;
5636 weightderiv(thisMeasurementIdx, bremNoBase + bremno) = qoverpBremTerm;
5637 for (
int bremno2 = nbremupstream; bremno2 < bremno; bremno2++) {
5638 weightderiv(thisMeasurementIdx, bremNoBase + bremno2) = qoverpBremTerm - qoverpTerm;
5657 const int nMeas = (
int)
res.size();
5662 for (
int i = 0;
i < nPerPars;
i++) {
5671 const std::unique_ptr<GXFTrackState> & state = trajectory.
trackStates()[
i];
5674 if (meff ==
nullptr) {
5675 measno += state->numberOfMeasuredParameters();
5679 const int firstMeasurement =
i < nUpstreamStates ? 0 : measno;
5680 const int lastMeasurement =
i < nUpstreamStates ? measno : nMeas - nBrem;
5683 && (trajectory.
prefit() == 0 || meff->
deltaE() == 0)) {
5684 const int scatterPos = nPerPars + 2 * scatno;
5696 const int bremPos = nPerPars + nScatPars + bremno;
5707 const Cache & cache,
5720 const int nMeas = (
int)
res.size();
5722 for (
int k = 0;
k < nFitPars;
k++) {
5731 for (
int measno = minMeasK; measno < maxMeasK; measno++) {
5732 b[
k] +=
res[measno] * (1. /
error[measno]) * weightDeriv(measno,
k);
5742 if (
k == 4 ||
k >= nPerPars + nScatPars) {
5743 for (
int measno = nMeas - nBrem; measno < nMeas; measno++) {
5744 b[
k] +=
res[measno] * (1. /
error[measno]) * weightDeriv(measno,
k);
5751 const Cache & cache,
5758 for (
int k = 0;
k < nFitPars;
k++) {
5759 for (
int l =
k;
l < nFitPars;
l++) {
5764 for (
int measno = minMeas; measno < maxMeas; measno++) {
5765 a_kl += weightDeriv(measno,
k) * weightDeriv(measno,
l);
5768 a.fillSymmetric(
l,
k, a_kl);
5786 const int nMeas = (
int)
res.size();
5793 for (
int k = nPerPars;
k < nPerPars + nScatPars;
k += 2) {
5803 for (
int measno = nMeas - nBrem; measno < nMeas; measno++) {
5804 for (
int k = 4;
k < nFitPars;
k++) {
5806 k = nPerPars + nScatPars;
5809 for (
int l =
k;
l < nFitPars;
l++) {
5811 l = nPerPars + nScatPars;
5814 const double a_kl =
a(
l,
k) + weightDeriv(measno,
k) * weightDeriv(measno,
l);
5815 a.fillSymmetric(
l,
k, a_kl);
5827 const double oldRedChi2,
5828 const double newRedChi2
5836 bool weightChanged =
false;
5843 double newPhiWeight = 1.1;
5844 double newThetaWeight = 1.001;
5845 if (trajectory.
prefit() == 0) {
5851 newPhiWeight = 1.00000001;
5852 }
else if (
it == 1) {
5853 newPhiWeight = 1.0000001;
5854 }
else if (
it <= 3) {
5855 newPhiWeight = 1.0001;
5856 }
else if (
it <= 6) {
5857 newPhiWeight = 1.01;
5860 if (newRedChi2 > oldRedChi2 - 1 && newRedChi2 < oldRedChi2) {
5861 newPhiWeight = 1.0001;
5862 newThetaWeight = 1.0001;
5863 }
else if (newRedChi2 > oldRedChi2 - 25 && newRedChi2 < oldRedChi2) {
5864 newPhiWeight = 1.001;
5865 newThetaWeight = 1.0001;
5872 std::size_t scatno = 0;
5877 for (
const auto & state : trajectory.
trackStates()) {
5880 if (meff ==
nullptr) {
5884 const bool isValidPlaneSurface =
5886 static_cast<const PlaneSurface *
>(&state->associatedSurface()) !=
nullptr;
5891 if (meff->
deltaE() == 0 || (trajectory.
prefit() == 0 && isValidPlaneSurface)) {
5892 weightChanged =
true;
5894 const int scatNoIndex = 2 * scatno + nPerPars;
5899 message <<
"scatno is out of range " << scatno <<
" !< " << cache.
m_phiweight.size();
5900 throw std::range_error(
message.str());
5908 a(scatNoIndex, scatNoIndex) /= cache.
m_phiweight[scatno];
5912 a(scatNoIndex, scatNoIndex) *= newPhiWeight;
5913 }
else if (trajectory.
prefit() >= 2) {
5914 a(scatNoIndex, scatNoIndex) *= newPhiWeight;
5915 a(scatNoIndex + 1, scatNoIndex + 1) *= newThetaWeight;
5942 trajectory.
prefit() == 2 &&
5945 (newRedChi2 < oldRedChi2 - 25 || newRedChi2 > oldRedChi2)
5950 return weightChanged;
5959 std::size_t scatno = 0;
5970 message <<
"scatno is out of range " << scatno <<
" !< " << cache.
m_phiweight.size();
5971 throw std::range_error(
message.str());
5974 const bool isValidPlaneSurface =
5976 static_cast<const PlaneSurface *
>(&state->associatedSurface()) !=
nullptr;
5978 if (meff->
deltaE() == 0 || isValidPlaneSurface) {
5979 const int scatNoIndex = 2 * scatno + nPerPars;
5980 a(scatNoIndex, scatNoIndex) /= cache.
m_phiweight[scatno];
5995 const EventContext& ctx,
6004 const int nDOFold = trajectory.
nDOF();
6005 const double oldChi2 = trajectory.
chi2();
6006 const double oldRedChi2 = nDOFold > 0 ? oldChi2 / nDOFold : 0;
6028 int bremno_maxbrempull = 0;
6044 if ((state_maxbrempull !=
nullptr) && trajectory.
converged()) {
6053 const int nDOFnew = trajectory.
nDOF();
6054 const double newChi2 = trajectory.
chi2();
6055 const double newRedChi2 = nDOFnew > 0 ? newChi2 / nDOFnew : 0;
6057 ATH_MSG_DEBUG(
"old chi2: " << oldChi2 <<
"/" << nDOFold <<
"=" << oldRedChi2 <<
6058 ", new chi2: " << newChi2 <<
"/" << nDOFnew <<
"=" << newRedChi2);
6093 if (doDeriv || weightChanged) {
6104 if (trajectory.
prefit() == 0) {
6110 if (
nSiHits + nTrtHits != nHits) {
6117 (newRedChi2 < 2 || (newRedChi2 < oldRedChi2 && newRedChi2 > oldRedChi2 - .5))
6139 Eigen::LLT<Eigen::MatrixXd>
const llt(lu_m);
6141 if (llt.info() != Eigen::Success) {
6165 double d0 = refpar->parameters()[
Trk::d0];
6166 double z0 = refpar->parameters()[
Trk::z0];
6169 double qoverp = refpar->parameters()[
Trk::qOverP];
6171 if (nperparams > 0) {
6172 d0 += deltaParameters[0];
6173 z0 += deltaParameters[1];
6174 phi += deltaParameters[2];
6175 theta += deltaParameters[3];
6176 qoverp = (trajectory.
m_straightline) ? 0 : .001 * deltaParameters[4] + qoverp;
6188 std::vector < std::pair < double, double >>&scatangles = trajectory.
scatteringAngles();
6189 for (
int i = 0;
i < nscat;
i++) {
6190 scatangles[
i].first += deltaParameters[2 *
i + nperparams];
6191 scatangles[
i].second += deltaParameters[2 *
i + nperparams + 1];
6197 std::vector < double >&delta_ps = trajectory.
brems();
6198 for (
int i = 0;
i < nbrem;
i++) {
6199 delta_ps[
i] += deltaParameters[nperparams + 2 * nscat +
i];
6205 std::unique_ptr<const TrackParameters> newper(
6225 const EventContext& evtctx
6236 if (!splitProbContainer.
isValid()) {
6240 std::vector<std::unique_ptr<GXFTrackState>> &
states = trajectory.
trackStates();
6247 for (
size_t stateno = 0; stateno <
states.size(); stateno++) {
6252 measno +=
states[stateno-1]->numberOfMeasuredParameters();
6255 std::unique_ptr<GXFTrackState> & state =
states[stateno];
6267 const auto *
const rot =
static_cast<const RIO_OnTrack *
>(pMeas);
6279 if (!splitProb.isSplit()) {
6280 ATH_MSG_DEBUG(
"Pixel cluster is not split so no need to update" );
6284 std::unique_ptr < const RIO_OnTrack > newrot;
6285 double *olderror = state->measurementErrors();
6288 double newerror[5] = {-1,-1,-1,-1,-1};
6289 double newres[2] = {-1,-1};
6291 newrot.reset(
m_ROTcreator->correct(*prd, *trackpars, evtctx));
6298 newerror[0] = std::sqrt(
covmat(0, 0));
6299 newres[0] = newrot->localParameters()[
Trk::locX] - trackpars->parameters()[
Trk::locX];
6300 newerror[1] = std::sqrt(
covmat(1, 1));
6301 newres[1] = newrot->localParameters()[
Trk::locY] - trackpars->parameters()[
Trk::locY];
6303 if (
a.cols() != nfitpars) {
6308 for(
int k =0;
k<2;
k++ ){
6309 const double oldres =
res[measno+
k];
6310 res[measno+
k] = newres[
k];
6311 err[measno+
k] = newerror[
k];
6313 for (
int i = 0;
i < nfitpars;
i++) {
6314 if (weightderiv(measno+
k,
i) == 0) {
6318 b[
i] -= weightderiv(measno+
k,
i) * (oldres / olderror[
k] - (newres[
k] * olderror[
k]) / (newerror[
k] * newerror[
k]));
6320 for (
int j =
i; j < nfitpars; j++) {
6324 weightderiv(measno+
k,
i) *
6325 weightderiv(measno+
k, j) *
6326 ((olderror[
k] * olderror[
k]) / (newerror[
k] * newerror[
k]) - 1)
6330 weightderiv(measno+
k,
i) *= olderror[
k] / newerror[
k];
6334 state->setMeasurement(std::move(newrot));
6335 state->setMeasurementErrors(newerror);
6350 const EventContext& ctx
6358 std::vector<std::unique_ptr<GXFTrackState>> &
states = trajectory.
trackStates();
6364 if (
a.cols() != nfitpars) {
6372 bool outlierremoved =
false;
6373 bool hitrecalibrated =
false;
6375 for (
int stateno = 0; stateno < (
int)
states.size(); stateno++) {
6376 std::unique_ptr<GXFTrackState> & state =
states[stateno];
6384 std::abs(state->trackParameters()->parameters()[
Trk::driftRadius]) > 1.05 * state->associatedSurface().bounds().r()
6389 outlierremoved =
true;
6391 double *
errors = state->measurementErrors();
6392 const double olderror =
errors[0];
6396 for (
int i = 0;
i < nfitpars;
i++) {
6397 if (weightderiv(measno,
i) == 0) {
6401 b[
i] -=
res[measno] * weightderiv(measno,
i) / olderror;
6403 for (
int j =
i; j < nfitpars; j++) {
6406 a(
i, j) - weightderiv(measno,
i) * weightderiv(measno, j)
6409 weightderiv(measno,
i) = 0;
6413 }
else if (trtrecal) {
6414 double *
errors = state->measurementErrors();
6415 const double olderror =
errors[0];
6417 const auto *
const thisMeasurement{state->measurement()};
6423 if (oldrot->prepRawData() !=
nullptr) {
6424 const double dcradius = oldrot->prepRawData()->localPosition()[
Trk::driftRadius];
6426 const double trackradius = state->trackParameters()->parameters()[
Trk::driftRadius];
6428 std::unique_ptr<const Trk::RIO_OnTrack> newrot =
nullptr;
6429 const double distance = std::abs(std::abs(trackradius) - dcradius);
6431 if (
distance < scalefactor * dcerror && (olderror > 1. || trackradius * oldradius < 0)) {
6432 newrot.reset(
m_ROTcreator->correct(*oldrot->prepRawData(), *state->trackParameters(), ctx));
6433 }
else if (
distance > scalefactor * dcerror && olderror < 1.) {
6434 newrot.reset(
m_broadROTcreator->correct(*oldrot->prepRawData(), *state->trackParameters(), ctx));
6437 if (newrot !=
nullptr) {
6439 hitrecalibrated =
true;
6443 if ((measno < 0) or (measno >= (
int)
res.size())) {
6444 throw std::runtime_error(
6445 "'res' array index out of range in TrkGlobalChi2Fitter/src/GlobalChi2Fitter.cxx:" +
std::to_string(__LINE__)
6449 const double oldres =
res[measno];
6450 const double newres = newradius - state->trackParameters()->parameters()[
Trk::driftRadius];
6452 state->setMeasurement(std::move(newrot));
6456 for (
int i = 0;
i < nfitpars;
i++) {
6457 if (weightderiv(measno,
i) == 0) {
6461 b[
i] -= weightderiv(measno,
i) * (oldres / olderror - (newres * olderror) / (newerror * newerror));
6463 for (
int j =
i; j < nfitpars; j++) {
6470 i < nperpars + 2 * nscats &&
6471 (
i - nperpars) % 2 == 0
6478 a(
i, j) + weightderiv(measno,
i) * weightderiv(measno, j) * ((olderror * olderror) / (newerror * newerror) - 1) *
weight
6481 weightderiv(measno,
i) *= olderror / newerror;
6484 res[measno] = newres;
6485 err[measno] = newerror;
6494 measno += state->numberOfMeasuredParameters();
6498 if (trajectory.
nDOF() < 0) {
6503 if (outlierremoved || hitrecalibrated) {
6512 const EventContext& ctx,
6520 bool trackok =
false;
6522 std::unique_ptr < GXFTrajectory > cleanup_oldtrajectory;
6524 std::unique_ptr < GXFTrajectory > cleanup_newtrajectory;
6530 while (!trackok && oldtrajectory->
nDOF() > 0) {
6532 std::vector<std::unique_ptr<GXFTrackState>> &
states = oldtrajectory->
trackStates();
6540 if (nhits != nsihits) {
6544 double maxsipull = -1;
6546 int hitno_maxsipull = -1;
6547 int measno_maxsipull = -1;
6548 int stateno_maxsipull = 0;
6560 for (
int stateno = 0; stateno < (
int)
states.size(); stateno++) {
6561 std::unique_ptr<GXFTrackState> & state =
states[stateno];
6567 double *
errors = state->measurementErrors();
6569 const Amg::MatrixX & hitcov = state->measurement()->localCovariance();
6570 const double sinstereo = state->sinStereo();
6571 const double cosstereo = (sinstereo == 0) ? 1 : std::sqrt(1 - sinstereo * sinstereo);
6572 double weight1 = -1;
6574 if (hitcov(0, 0) > trackcov(0, 0)) {
6575 if (sinstereo == 0) {
6579 trackcov(0, 0) * cosstereo * cosstereo + 2 *
6580 trackcov(1, 0) * cosstereo * sinstereo + trackcov(1, 1) * sinstereo * sinstereo
6585 const double weight2 = (
6591 double sipull1 = weight1 > 0 ? std::abs(
res[measno] / std::sqrt(weight1)) : -1;
6592 const double sipull2 = (
6594 std::abs(
res[measno + 1] / std::sqrt(weight2)) :
6597 sipull1 =
std::max(sipull1, sipull2);
6599 if (sipull1 > maxsipull) {
6600 maxsipull = sipull1;
6601 measno_maxsipull = measno;
6602 state_maxsipull = state.get();
6603 stateno_maxsipull = stateno;
6604 hitno_maxsipull = hitno;
6615 measno += state->numberOfMeasuredParameters();
6619 const double maxpull = maxsipull;
6621 ATH_MSG_DEBUG(
" maxsipull: " << maxsipull <<
" hitno_maxsipull: " <<
6622 hitno_maxsipull <<
" n3sigma: " << n3sigma <<
" cut: " <<
cut <<
" cut2: " << cut2);
6633 state_maxsipull = oldtrajectory->
trackStates()[stateno_maxsipull].get();
6636 const auto *
const rot =
static_cast<const RIO_OnTrack *
>(pMeas);
6639 std::unique_ptr < const RIO_OnTrack > broadrot;
6644 Amg::VectorX parameterVector = trackpar_maxsipull->parameters();
6645 const std::unique_ptr<const TrackParameters> trackparForCorrect(
6654 state_maxsipull->trackCovariance())
6658 newerror[0] = newerror[1] = newerror[2] = newerror[3] = newerror[4] = -1;
6659 double newpull = -1;
6660 double newpull1 = -1;
6661 double newpull2 = -1;
6662 double newres1 = -1;
6663 double newres2 = -1;
6664 double newsinstereo = 0;
6679 if (state_maxsipull->
sinStereo() != 0) {
6680 const auto [covEigenValueSmall, covStereoAngle] = principalComponentAnalysis2x2(
covmat);
6681 newerror[0] = std::sqrt(covEigenValueSmall);
6682 newsinstereo =
std::sin(covStereoAngle);
6684 newerror[0] = std::sqrt(
covmat(0, 0));
6687 const double cosstereo = (newsinstereo == 0) ? 1. : std::sqrt(1 - newsinstereo * newsinstereo);
6689 if (cosstereo != 1.) {
6691 cosstereo * (broadrot->localParameters()[
Trk::locX] - trackpar_maxsipull->parameters()[
Trk::locX]) +
6692 newsinstereo * (broadrot->localParameters()[
Trk::locY] - trackpar_maxsipull->parameters()[
Trk::locY])
6695 newres1 = broadrot->localParameters()[
Trk::locX] - trackpar_maxsipull->parameters()[
Trk::locX];
6698 if (newerror[0] == 0.0) {
6699 ATH_MSG_WARNING(
"Measurement error is zero or negative, treating as outlier");
6702 newpull1 = std::abs(newres1 / newerror[0]);
6706 newerror[1] = std::sqrt(
covmat(1, 1));
6707 newres2 = broadrot->localParameters()[
Trk::locY] - trackpar_maxsipull->parameters()[
Trk::locY];
6708 newpull2 = std::abs(newres2 / newerror[1]);
6711 newpull =
std::max(newpull1, newpull2);
6717 (newerror[0] > 1.5 * olderror[0] || newerror[1] > 1.5 * std::abs(olderror[1]))
6719 if ((measno_maxsipull < 0) or(measno_maxsipull >= (
int)
res.size())) {
6720 throw std::runtime_error(
6721 "'res' array index out of range in TrkGlobalChi2Fitter/src/GlobalChi2Fitter.cxx:" +
std::to_string(__LINE__)
6726 newtrajectory = oldtrajectory;
6728 if (
a.cols() != nfitpars) {
6732 const double oldres1 =
res[measno_maxsipull];
6733 res[measno_maxsipull] = newres1;
6734 err[measno_maxsipull] = newerror[0];
6736 for (
int i = 0;
i < nfitpars;
i++) {
6737 if (weightderiv(measno_maxsipull,
i) == 0) {
6741 b[
i] -= weightderiv(measno_maxsipull,
i) * (oldres1 / olderror[0] - (newres1 * olderror[0]) / (newerror[0] * newerror[0]));
6743 for (
int j =
i; j < nfitpars; j++) {
6747 weightderiv(measno_maxsipull,
i) *
6748 weightderiv(measno_maxsipull, j) *
6749 ((olderror[0] * olderror[0]) / (newerror[0] * newerror[0]) - 1)
6753 weightderiv(measno_maxsipull,
i) *= olderror[0] / newerror[0];
6757 const double oldres2 =
res[measno_maxsipull + 1];
6758 res[measno_maxsipull + 1] = newres2;
6759 err[measno_maxsipull + 1] = newerror[1];
6761 for (
int i = 0;
i < nfitpars;
i++) {
6762 if (weightderiv(measno_maxsipull + 1,
i) == 0) {
6766 b[
i] -= weightderiv(measno_maxsipull + 1,
i) * (oldres2 / olderror[1] - (newres2 * olderror[1]) / (newerror[1] * newerror[1]));
6768 for (
int j =
i; j < nfitpars; j++) {
6772 weightderiv(measno_maxsipull + 1,
i) *
6773 weightderiv(measno_maxsipull + 1, j) *
6774 ((olderror[1] * olderror[1]) / (newerror[1] * newerror[1]) - 1)
6779 weightderiv(measno_maxsipull + 1,
i) *= olderror[1] / newerror[1];
6784 "Recovering outlier, hitno=" << hitno_maxsipull <<
" measno=" <<
6785 measno_maxsipull <<
" pull=" << maxsipull <<
" olderror_0=" <<
6786 olderror[0] <<
" newerror_0=" << newerror[0] <<
" olderror_1=" <<
6787 olderror[1] <<
" newerror_1=" << newerror[1]
6796 ((n3sigma < 2 && maxsipull > cut2 && maxsipull <
cut) || n3sigma > 1) &&
6806 "Removing outlier, hitno=" << hitno_maxsipull <<
", measno=" <<
6807 measno_maxsipull <<
" pull=" << maxsipull
6814 cleanup_newtrajectory = std::make_unique<GXFTrajectory>(*oldtrajectory);
6815 newtrajectory = cleanup_newtrajectory.get();
6817 if (newa.cols() != nfitpars) {
6823 if ((measno_maxsipull < 0) or(measno_maxsipull >= (
int)
res.size())) {
6824 throw std::runtime_error(
6825 "'res' array index out of range in TrkGlobalChi2Fitter/src/GlobalChi2Fitter.cxx:" +
std::to_string(__LINE__)
6829 const double oldres1 =
res[measno_maxsipull];
6830 newres[measno_maxsipull] = 0;
6832 for (
int i = 0;
i < nfitpars;
i++) {
6833 if (weightderiv(measno_maxsipull,
i) == 0) {
6837 newb[
i] -= weightderiv(measno_maxsipull,
i) * oldres1 / olderror[0];
6839 for (
int j =
i; j < nfitpars; j++) {
6843 weightderiv(measno_maxsipull,
i) *
6844 weightderiv(measno_maxsipull, j)
6848 newweightderiv(measno_maxsipull,
i) = 0;
6852 const double oldres2 =
res[measno_maxsipull + 1];
6853 newres[measno_maxsipull + 1] = 0;
6855 for (
int i = 0;
i < nfitpars;
i++) {
6856 if (weightderiv(measno_maxsipull + 1,
i) == 0) {
6860 newb[
i] -= weightderiv(measno_maxsipull + 1,
i) * oldres2 / olderror[1];
6862 for (
int j =
i; j < nfitpars; j++) {
6863 if (weightderiv(measno_maxsipull + 1, j) == 0) {
6870 weightderiv(measno_maxsipull + 1,
i) *
6871 weightderiv(measno_maxsipull + 1, j)
6875 newweightderiv(measno_maxsipull + 1,
i) = 0;
6879 newtrajectory->
setOutlier(stateno_maxsipull);
6903 ctx, cache, *newtrajectory,
it, *newap, *newbp, lu_m, doderiv);
6919 const double oldchi2 = oldtrajectory->
chi2() / oldtrajectory->
nDOF();
6920 const double newchi2 = (newtrajectory->
nDOF() > 0) ? newtrajectory->
chi2() / newtrajectory->
nDOF() : 0;
6923 if (newtrajectory->
nDOF() != oldtrajectory->
nDOF() && maxsipull > cut2) {
6924 mindiff = (oldchi2 > .33 *
m_chi2cut || noutl > 0) ? .8 : 1.;
6926 if (noutl == 0 && maxsipull <
cut - .5 && oldchi2 < .5 *
m_chi2cut) {
6931 if (newchi2 > oldchi2 || (newchi2 > oldchi2 - mindiff && newchi2 > .33 * oldchi2)) {
6932 ATH_MSG_DEBUG(
"Outlier not confirmed, keeping old trajectory");
6940 (void)cleanup_oldtrajectory.release();
6941 return oldtrajectory;
6943 if (oldtrajectory != newtrajectory) {
6944 cleanup_oldtrajectory = std::move(cleanup_newtrajectory);
6945 oldtrajectory = newtrajectory;
6952 Eigen::LLT < Eigen::MatrixXd >
const lltOfW(
a);
6953 if (lltOfW.info() == Eigen::Success) {
6957 const int ncols =
a.cols();
6958 Amg::MatrixX const weightInvAMG = Amg::MatrixX::Identity(ncols, ncols);
6959 fullcov = lltOfW.solve(weightInvAMG);
6977 oldtrajectory->
nDOF() > 0 &&
6986 (void)cleanup_oldtrajectory.release();
6987 return oldtrajectory;
6999 if (
const auto *pMeas{hit->measurement()};
7005 nrealmeas += hit->numberOfMeasuredParameters();
7015 if (
const auto *pMeas{hit->measurement()};
7021 for (
int i = measindex;
i < measindex + hit->numberOfMeasuredParameters();
i++) {
7024 if ((j == 4 && !oldtrajectory.
m_straightline) || j >= nperpars + 2 * nscat) {
7032 measindex += hit->numberOfMeasuredParameters();
7033 }
else if (hit->materialEffects() ==
nullptr) {
7034 measindex2 += hit->numberOfMeasuredParameters();
7040 const EventContext & ctx,
7048 std::unique_ptr<const TrackParameters> per(
nullptr);
7051 std::unique_ptr<const TrackParameters> prevpar(
7056 std::vector<std::pair<const Layer *, const Layer *>> & upstreamlayers = oldtrajectory.
upstreamMaterialLayers();
7060 if (prevpar ==
nullptr) {
7065 const Layer *
layer = layer1 !=
nullptr ? layer1 : layer2;
7068 prevpar->position(), prevpar->momentum().unit()
7086 std::unique_ptr<const TrackParameters> layerpar(
7090 layer->surfaceRepresentation(),
7098 if (layerpar ==
nullptr) {
7102 if (
layer->surfaceRepresentation().bounds().inside(layerpar->localPosition())) {
7106 prevpar = std::move(layerpar);
7119 if (startfactor > 0.5) {
7120 std::unique_ptr<const TrackParameters> updatedpar =
m_matupdator->update(
7124 if (updatedpar !=
nullptr) {
7144 if (endfactor > 0.5) {
7145 std::unique_ptr<const TrackParameters> updatedpar =
m_matupdator->update(
7149 if (updatedpar !=
nullptr) {
7155 if (prevpar !=
nullptr) {
7167 if (per ==
nullptr) {
7168 ATH_MSG_DEBUG(
"Failed to extrapolate to perigee, returning 0");
7187 std::unique_ptr<GXFTrackState>
7189 const EventContext & ctx,
7196 if (per ==
nullptr) {
7200 ATH_MSG_DEBUG(
"Final perigee: " << *per <<
" pos: " << per->position() <<
" pT: " << per->pT());
7206 const std::vector<std::unique_ptr<TrackParameters>> & hc,
7207 std::set<Identifier> & id_set,
7208 std::set<Identifier> & sct_set,
7218 for (
const std::unique_ptr<TrackParameters> &
tp : hc) {
7224 if (
tp ==
nullptr) {
7235 if (de ==
nullptr) {
7246 if (id_set.find(
id) != id_set.end()) {
7295 if (sct_set.find(
os) != sct_set.end()) {
7296 ++
rv.m_sct_double_hole;
7325 for (
const std::unique_ptr<GXFTrackState> &
s : trajectory.
trackStates()) {
7337 std::vector<std::reference_wrapper<GXFTrackState>>
rv;
7345 for (
const std::unique_ptr<GXFTrackState> &
s : trajectory.
trackStates()) {
7359 rv.emplace_back(*
s);
7366 const TrkDetElementBase * de =
s->trackParameters()->associatedSurface().associatedDetectorElement();
7368 if (de !=
nullptr) {
7381 if (
s.get() == lastmeas) {
7391 const EventContext & ctx,
7392 const std::vector<std::reference_wrapper<GXFTrackState>> &
states
7404 constexpr
uint min_meas = 3;
7405 if (std::count_if(
states.begin(),
states.end(), [](
const GXFTrackState &
s){ return s.getStateType(TrackStateOnSurface::Measurement); }) < min_meas) {
7409 bool seen_meas =
false;
7411 std::set<Identifier> id_set;
7412 std::set<Identifier> sct_set;
7418 for (std::size_t
i = 0;
i <
states.size() - 1;
i++) {
7441 const double dist = (
beg.trackParameters()->position() -
end.trackParameters()->position()).
norm();
7443 const bool zStartValid = std::abs(
beg.trackParameters()->position().z())<10000.;
7445 ATH_MSG_DEBUG(
"Pathological track parameter well outside of detector");
7446 ATH_MSG_DEBUG(
"Propagator might have issue with this, skipping");
7456 if (seen_meas && dist >= 2.5 && zStartValid) {
7463 std::optional<std::vector<std::unique_ptr<TrackParameters>>> & hc =
beg.getHoles();
7464 std::vector<std::unique_ptr<TrackParameters>>
states;
7472 if (hc.has_value()) {
7511 std::vector<std::unique_ptr<Trk::TrackParameters>>
const bl =
m_extrapolator->extrapolateBlindly(
7532 const EventContext & ctx,
7538 auto trajectory = std::make_unique<Trk::TrackStates>();
7546 std::unique_ptr<GXFTrackState> perigee_ts =
makeTrackFindPerigee(ctx, cache, oldtrajectory, matEffects);
7548 if (perigee_ts ==
nullptr) {
7554 trajectory->reserve(tmptrajectory.
trackStates().size());
7560 hit->resetTrackCovariance();
7566 hit->materialEffects()))
7572 auto trackState = hit->trackStateOnSurface();
7573 hit->resetTrackCovariance();
7574 trajectory->emplace_back(trackState.release());
7577 auto qual = std::make_unique<FitQuality>(tmptrajectory.
chi2(), tmptrajectory.
nDOF());
7597 std::unique_ptr<Track>
rv = std::make_unique<Track>(
info, std::move(trajectory), std::move(
qual));
7607 std::unique_ptr<TrackSummary>
ts = std::make_unique<TrackSummary>();
7616 std::optional<TrackHoleCount> hole_count;
7641 if (hole_count.has_value()) {
7654 rv->setTrackSummary(std::move(
ts));
7660 GlobalChi2Fitter::~GlobalChi2Fitter() =
default;
7663 const EventContext & ctx,
7672 std::vector<std::unique_ptr<TrackParameters>>
rv =
m_extrapolator->extrapolateStepwise(
7686 &
rv.front()->associatedSurface() == &
src.associatedSurface() ||
7687 trackParametersClose(*
rv.front(),
src, 0.001) ||
7691 rv.front().reset(
nullptr);
7702 &
rv.back()->associatedSurface() == &
src.associatedSurface() ||
7703 trackParametersClose(*
rv.back(),
src, 0.001) ||
7707 rv.back().reset(
nullptr);
7714 const EventContext & ctx,
7722 std::unique_ptr<const TrackParameters>
rv;
7723 std::optional<TransportJacobian> jac{};
7734 if (
rv !=
nullptr && calcderiv) {
7739 std::optional<std::vector<std::unique_ptr<TrackParameters>>> extrapolation;
7748 std::move(extrapolation)
7753 const EventContext & ctx,
7764 ctx, prev,
ts, propdir, bf, calcderiv, holesearch
7767 if (
rv.m_parameters ==
nullptr) {
7768 propdir = invertPropdir(propdir);
7771 ctx, prev,
ts, propdir, bf, calcderiv, holesearch
7779 const EventContext& ctx,
7786 std::vector<std::unique_ptr<GXFTrackState>> &
states = trajectory.
trackStates();
7789 std::unique_ptr<const TrackParameters> tmptrackpar;
7791 for (
int hitno = nstatesupstream - 1; hitno >= 0; hitno--) {
7821 (
rv.m_parameters !=
nullptr) &&
7822 (prevtrackpar->
position() -
rv.m_parameters->position()).mag() > 5 *
mm
7828 if (
rv.m_parameters ==
nullptr) {
7829 ATH_MSG_DEBUG(
"propagation failed, prev par: " << *prevtrackpar <<
7830 " pos: " << prevtrackpar->
position() <<
" destination surface: " << surf1);
7834 states[hitno]->setTrackParameters(std::move(
rv.m_parameters));
7838 if (
rv.m_jacobian != std::nullopt) {
7840 states[hitno]->materialEffects() !=
nullptr &&
7841 states[hitno]->materialEffects()->deltaE() != 0 &&
7842 states[hitno]->materialEffects()->sigmaDeltaE() <= 0 &&
7845 const double p = 1. / std::abs(currenttrackpar->parameters()[
Trk::qOverP]);
7846 const double de = std::abs(
states[hitno]->materialEffects()->deltaE());
7847 const double mass = trajectory.
mass();
7848 const double newp = std::sqrt(
p *
p + 2 * de * std::sqrt(
mass *
mass +
p *
p) + de * de);
7849 (*
rv.m_jacobian) (4, 4) = ((
p +
p * de / std::sqrt(
p *
p +
mass *
mass)) / newp) *
p *
p / (newp * newp);
7852 states[hitno]->setJacobian(*
rv.m_jacobian);
7853 }
else if (calcderiv) {
7860 if (meff !=
nullptr && hitno != 0) {
7862 surf, *meff, *
states[hitno]->trackParameters(), trajectory.
mass(), -1
7865 if (std::holds_alternative<FitterStatusCode>(
r)) {
7866 return std::get<FitterStatusCode>(
r);
7869 tmptrackpar = std::move(
std::get<std::unique_ptr<const TrackParameters>>(
r));
7870 prevtrackpar = tmptrackpar.get();
7872 prevtrackpar = currenttrackpar;
7878 for (
int hitno = nstatesupstream; hitno < (
int)
states.size(); hitno++) {
7900 (
rv.m_parameters !=
nullptr) &&
7902 (prevtrackpar->
position() -
rv.m_parameters->position()).mag() > 5 *
mm
7907 if (
rv.m_parameters ==
nullptr) {
7908 ATH_MSG_DEBUG(
"propagation failed, prev par: " << *prevtrackpar <<
7909 " pos: " << prevtrackpar->
7910 position() <<
" destination surface: " << surf);
7914 if (
rv.m_jacobian != std::nullopt) {
7916 states[hitno]->materialEffects() !=
nullptr &&
7917 states[hitno]->materialEffects()->deltaE() != 0 &&
7918 states[hitno]->materialEffects()->sigmaDeltaE() <= 0 &&
7921 const double p = 1 / std::abs(
rv.m_parameters->parameters()[
Trk::qOverP]);
7922 const double de = std::abs(
states[hitno]->materialEffects()->deltaE());
7923 const double mass = trajectory.
mass();
7924 double newp =
p *
p - 2 * de * std::sqrt(
mass *
mass +
p *
p) + de * de;
7927 newp = std::sqrt(newp);
7930 (*
rv.m_jacobian) (4, 4) = ((
p -
p * de / std::sqrt(
p *
p +
mass *
mass)) / newp) *
p *
p / (newp * newp);
7933 states[hitno]->setJacobian(*
rv.m_jacobian);
7934 }
else if (calcderiv) {
7941 if (meff !=
nullptr) {
7943 surf, *meff, *
rv.m_parameters, trajectory.
mass(), +1
7946 if (std::holds_alternative<FitterStatusCode>(
r)) {
7947 return std::get<FitterStatusCode>(
r);
7950 rv.m_parameters = std::move(
std::get<std::unique_ptr<const TrackParameters>>(
r));
7953 states[hitno]->setTrackParameters(std::move(
rv.m_parameters));
7954 prevtrackpar =
states[hitno]->trackParameters();
7973 ATH_MSG_DEBUG(
"Angles out of range, phi: " << newphi <<
" theta: " << newtheta);
7977 double newqoverp = 0;
7991 newqoverp = std::copysign(1 / std::sqrt(newp2),
old[
Trk::qOverP]);
7998 old[0],
old[1], newphi, newtheta, newqoverp, std::nullopt
8010 using Matrix55 = Eigen::Matrix<double, 5, 5>;
8012 Matrix55 initialjac;
8013 initialjac.setZero();
8014 initialjac(4, 4) = 1;
8016 Matrix55 jacvertex(initialjac);
8018 std::vector<Matrix55, Eigen::aligned_allocator<Matrix55>> jacscat(trajectory.
numberOfScatterers(), initialjac);
8019 std::vector<Matrix55, Eigen::aligned_allocator<Matrix55>> jacbrem(trajectory.
numberOfBrems(), initialjac);
8021 std::vector<std::unique_ptr<GXFTrackState>> &
states = trajectory.
trackStates();
8030 for (
const bool forward : {
false,
true}) {
8032 hit_begin = nstatesupstream;
8034 scatno = nscatupstream;
8035 bremno = nbremupstream;
8037 hit_begin = nstatesupstream - 1;
8044 int hitno = hit_begin;
8045 forward ? (hitno < hit_end) : (hitno >= hit_end);
8046 hitno += (forward ? 1 : -1)
8049 state =
states[hitno].get();
8053 if (fillderivmat && state->
derivatives().cols() != nfitpars) {
8061 const int jmaxbrem = 4;
8063 if (hitno == (forward ? hit_end - 1 : 0)) {
8064 if (!fillderivmat) {
8072 Eigen::Matrix<double, 5, 5> & jac = state->
jacobian();
8074 if (hitno == nstatesupstream + (forward ? 0 : -1)) {
8075 jacvertex.block<4, 5>(0, 0) = jac.block<4, 5>(0, 0);
8076 jacvertex(4, 4) = jac(4, 4);
8086 jcnt = jmax - jmin + 1;
8088 lp_bgn = forward ? nscatupstream : nscatupstream - 1;
8091 for (
int i = lp_bgn; forward ? (
i < lp_end) : (
i > lp_end);
i += (forward ? 1 : -1)) {
8093 i == scatno + (forward ? -1 : 1) &&
8094 prevstate !=
nullptr &&
8098 jacscat[
i].block(0, jmin, 4, jcnt) = jac.block(0, jmin, 4, jcnt);
8099 jacscat[
i](4, 4) = jac(4, 4);
8101 calculateJac(jac, jacscat[
i], jmin, jmax);
8105 Eigen::MatrixXd & derivmat = state->
derivatives();
8106 const int scatterPos = nperpars + 2 *
i;
8108 derivmat.block<4, 2>(0, scatterPos) = (forward ? 1 : -1) * jacscat[
i].block<4, 2>(0, 2);
8114 jcnt = jmax - jmin + 1;
8116 lp_bgn = forward ? nbremupstream : nbremupstream - 1;
8119 for (
int i = lp_bgn; forward ? (
i < lp_end) : (
i > lp_end);
i += (forward ? 1 : -1)) {
8121 i == bremno + (forward ? -1 : 1) &&
8126 jacbrem[
i].block(0, jmin, 4, jcnt) = jac.block(0, jmin, 4, jcnt);
8127 jacbrem[
i](4, 4) = jac(4, 4);
8129 calculateJac(jac, jacbrem[
i], jmin, jmax);
8133 Eigen::MatrixXd & derivmat = state->
derivatives();
8134 const int scatterPos = nperpars + 2 * nscats +
i;
8136 derivmat.block<5, 1>(0, scatterPos) = (forward ? .001 : -.001) * jacbrem[
i].block<5, 1>(0, 4);
8140 calculateJac(jac, jacvertex, 0, 4);
8144 Eigen::MatrixXd & derivmat = state->
derivatives();
8145 derivmat.block(0, 0, 4, nperpars) = jacvertex.block(0, 0, 4, nperpars);
8147 if (nperpars == 5) {
8148 derivmat.col(4).segment(0, 4) *= .001;
8149 derivmat(4, 4) = .001 * jacvertex(4, 4);
8155 (!trajectory.
prefit() ||
states[hitno]->materialEffects()->deltaE() == 0)
8157 scatno += (forward ? 1 : -1);
8161 states[hitno]->materialEffects() &&
8162 states[hitno]->materialEffects()->sigmaDeltaE() > 0
8164 bremno += (forward ? 1 : -1);
8167 prevstate =
states[hitno].get();
8176 bool onlylocal)
const {
8182 std::vector<std::unique_ptr<GXFTrackState>> &
states = trajectory.
trackStates();
8186 int i = nstatesupstream;
8187 for (
int j = 0; j < (
int)
states.size(); j++) {
8188 if (j < nstatesupstream) {
8195 for (
int stateno = 0; stateno < (
int)
states.size(); stateno++) {
8196 if (stateno == 0 || stateno == nstatesupstream) {
8197 prevstate =
nullptr;
8200 std::unique_ptr<GXFTrackState> & state =
states[
index];
8201 if (state->materialEffects() !=
nullptr) {
8202 prevstate = state.get();
8206 if (!state->hasTrackCovariance()) {
8207 state->zeroTrackCovariance();
8209 AmgMatrix(5, 5) & trackerrmat = state->trackCovariance();
8211 if ((prevstate !=
nullptr) &&
8215 Eigen::Matrix<double, 5, 5> & jac = state->jacobian();
8218 trackerrmat = jac * prevcov * jac.transpose();
8222 trackerrmat = derivatives * fullcovmat * derivatives.transpose();
8232 bool errorok =
true;
8233 for (
int i = 0;
i < 5;
i++) {
8236 && trackerrmat(
i,
i) > meascov(j, j)) {
8238 const double scale = std::sqrt(meascov(j, j) / trackerrmat(
i,
i));
8239 trackerrmat(
i,
i) = meascov(j, j);
8240 for (
int k = 0;
k < 5;
k++) {
8250 for (
int i = 0;
i < 5;
i++) {
8254 for (
int j = 0; j < 5; j++) {
8262 trackerrmat(4, 4) = 1
e-20;
8266 state->trackParameters();
8268 std::optional<
AmgMatrix(5, 5)> trkerrmat;
8270 if (state->hasTrackCovariance()) {
8271 trkerrmat = (state->trackCovariance());
8273 trkerrmat = std::nullopt;
8276 const AmgVector(5) & tpars = tmptrackpar->parameters();
8277 std::unique_ptr<const TrackParameters> trackpar(
8283 std::move(trkerrmat))
8285 state->setTrackParameters(std::move(trackpar));
8288 if (errorok && trajectory.
nDOF() > 0) {
8289 fitQual =
m_updator->fullStateFitQuality(
8290 *state->trackParameters(),
8298 state->setFitQuality(fitQual);
8300 prevstate = state.get();
8304 std::optional<TransportJacobian>
8306 const EventContext& ctx,
8319 std::optional<TransportJacobian> jac = std::make_optional<TransportJacobian>(J);
8322 0.01, 0.01, 0.00001, 0.00001, 0.000000001
8333 for (
int i = 0;
i < 5;
i++) {
8336 if (thisdiscsurf &&
i == 1) {
8342 if (
i == 0 && thiscylsurf) {
8344 }
else if (
i == 1 && thisdiscsurf) {
8350 std::unique_ptr<const TrackParameters> parpluseps(
8360 const std::unique_ptr<const TrackParameters> parminuseps(
8371 std::unique_ptr<const TrackParameters> newparpluseps(
8382 std::unique_ptr<const TrackParameters> newparminuseps(
8397 if (newparpluseps ==
nullptr) {
8409 if (newparminuseps ==
nullptr) {
8421 if ((newparpluseps ==
nullptr) || (newparminuseps ==
nullptr)) {
8425 for (
int j = 0; j < 5; j++) {
8429 if (j == 0 && cylsurf) {
8431 }
else if (j == 1 && discsurf) {
8435 (*jac) (j,
i) =
diff / (2 * eps[
i]);
8448 (
"Configure the minimum number of Iterations via jobOptions");
8470 auto nmeas1 = pDataVector->size();
8471 const auto *pLastValue = (*pDataVector)[nmeas1 - 1];
8479 testrot =
static_cast<const RIO_OnTrack *
>(pLastValue);
8481 if (lastMeasIsCompetingRIO){
8487 if (testrot ==
nullptr) {
8488 const auto *pPenultimate = (*pDataVector)[nmeas1 - 2];
8491 if(penultimateIsRIO){
8492 testrot =
static_cast<const RIO_OnTrack *
>(pPenultimate);
8494 if (penultimateIsCompetingRIO){
8503 (testrot !=
nullptr) &&
8520 if (cond_obj ==
nullptr) {
8529 std::stringstream
msg;
8531 throw std::runtime_error(
msg.str());