57 #include "CLHEP/Units/SystemOfUnits.h"
70 #include <Eigen/Dense>
71 #include <Eigen/StdVector>
82 return distsol.
first();
85 std::abs(distsol.
first()) < std::abs(distsol.
second()) ?
96 limitInversePValue(
double qOverP){
97 const double magnitude = std::abs(
qOverP);
99 constexpr
double maxP{100.*10e6*
MeV};
100 constexpr
double minP{1.e-3*
MeV};
101 constexpr
double lo {1./maxP};
102 constexpr
double hi {1./minP};
103 double limited = std::clamp(magnitude, lo, hi);
104 return std::copysign(limited,
qOverP);
108 std::pair<const Trk::TrackParameters *, const Trk::TrackParameters *> getFirstLastIdPar(
const Trk::Track &
track) {
114 while ((firstidpar ==
nullptr) && parit !=
track.trackParameters()->end()) {
116 ((**parit).covariance() !=
nullptr) &&
125 parit =
track.trackParameters()->end();
129 ((**parit).covariance() !=
nullptr) &&
134 }
while ((lastidpar ==
nullptr) && parit !=
track.trackParameters()->begin());
136 return std::make_pair(firstidpar, lastidpar);
151 std::abs(
a.parameters()[0] -
b.parameters()[0]) <
e &&
152 std::abs(
a.parameters()[1] -
b.parameters()[1]) <
e &&
153 std::abs(
a.parameters()[2] -
b.parameters()[2]) <
e
164 calculateJac(Eigen::Matrix<double, 5, 5> &jac,
165 Eigen::Matrix<double, 5, 5> &
out,
166 int jmin,
int jmax) {
170 out.block(0, 0, 4, jmin).setZero();
174 out.block(0, jmax + 1, 4, 5 - (jmax + 1)).setZero();
177 out(4, 4) = jac(4, 4);
187 std::pair<double, double> principalComponentAnalysis2x2(
const Amg::MatrixX &
mat) {
188 const double trace =
mat(0, 0) +
mat(1, 1);
189 const double diagonalProduct =
mat(0, 0) *
mat(1, 1);
190 const double mat01Sq =
mat(0, 1) *
mat(0, 1);
191 const double discriminant = std::sqrt(trace * trace - 4. * (diagonalProduct - mat01Sq));
193 const double eigenValueSmall = 0.5 * (trace -
discriminant);
194 const double stereoAngle = 0.5 * std::asin(2 *
mat(0, 1) / (-
discriminant));
196 return std::make_pair(eigenValueSmall, stereoAngle);
202 const std::string &
t,
203 const std::string &
n,
231 ATH_MSG_ERROR(
"Hole search requested but no boundary check tool provided.");
232 return StatusCode::FAILURE;
257 ATH_MSG_WARNING(
"FillDerivativeMatrix option selected, switching off acceleration!");
281 ATH_MSG_ERROR(
"Hole search requested but track summaries are disabled.");
282 return StatusCode::FAILURE;
287 return StatusCode::SUCCESS;
293 if (m_fit_status[
S_FITS] > 0) {
296 <<
" track fits failed because of a matrix inversion failure");
298 <<
" tracks were rejected by the outlier logic");
300 <<
" track fits failed because of a propagation failure");
302 <<
" track fits failed because of an invalid angle (theta/phi)");
304 <<
" track fits failed because the fit did not converge");
306 <<
" tracks did not pass the chi^2 cut");
308 <<
" tracks were killed by the energy loss update");
311 return StatusCode::SUCCESS;
316 std::unique_ptr<Track>
318 const EventContext& ctx,
324 ATH_MSG_DEBUG(
"--> entering GlobalChi2Fitter::fit(Track,Track,)");
339 const Track *indettrack = firstismuon ? &intrk2 : &intrk1;
340 const Track *muontrack = firstismuon ? &intrk1 : &intrk2;
342 bool measphi =
false;
349 rot = &crot->rioOnTrack(0);
359 double dotprod2 = measdir.dot(
362 if (std::abs(dotprod1) < .5 && std::abs(dotprod2) < .5) {
370 auto [firstidpar, lastidpar] = getFirstLastIdPar(*indettrack);
372 if ((firstidpar ==
nullptr) || (lastidpar ==
nullptr)) {
376 std::unique_ptr<const TrackParameters> parforcalo =
unique_clone(firstismuon ? firstidpar : lastidpar);
379 const AmgVector(5) & newpars = parforcalo->parameters();
381 parforcalo=parforcalo->associatedSurface().createUniqueTrackParameters(
382 newpars[0], newpars[1], newpars[2], newpars[3], 1 / 5000., std::nullopt
386 std::vector < MaterialEffectsOnTrack > calomeots;
389 calomeots =
m_calotool->extrapolationSurfacesAndEffects(
393 parforcalo->associatedSurface(),
406 if (calomeots.empty()) {
411 std::unique_ptr<Track> track;
421 double qoverpid = measperid !=
nullptr ? measperid->parameters()[
Trk::qOverP] : 0;
422 double qoverpmuon = measpermuon !=
nullptr ? measpermuon->parameters()[
Trk::qOverP] : 0;
424 const AmgSymMatrix(5) * errmatid = measperid !=
nullptr ? measperid->covariance() :
nullptr;
425 const AmgSymMatrix(5) * errmatmuon = measpermuon !=
nullptr ? measpermuon->covariance() :
nullptr;
429 (errmatid !=
nullptr) &&
430 (errmatmuon !=
nullptr) &&
435 double piderror = std::sqrt((*errmatid) (4, 4)) / (qoverpid * qoverpid);
436 double pmuonerror = std::sqrt((*errmatmuon) (4, 4)) / (qoverpmuon * qoverpmuon);
437 double energyerror = std::sqrt(
438 calomeots[1].energyLoss()->sigmaDeltaE() *
439 calomeots[1].energyLoss()->sigmaDeltaE() + piderror * piderror +
440 pmuonerror * pmuonerror
444 (std::abs(calomeots[1].energyLoss()->deltaE()) -
445 std::abs(1 / qoverpid) + std::abs(1 / qoverpmuon)
448 ATH_MSG_DEBUG(
"Changing from measured to parametrized energy loss");
453 parforcalo->associatedSurface(),
458 if (calomeots.empty()) {
466 bool firstfitwasattempted =
false;
469 if (!caloEntranceIsValid) {
480 qoverpid * qoverpmuon > 0
486 firstfitwasattempted =
true;
491 (track ==
nullptr) &&
492 !firstfitwasattempted &&
499 trajectory = trajectory2;
503 bool pseudoupdated =
false;
505 if (track !=
nullptr) {
506 for (std::unique_ptr<GXFTrackState> & pseudostate : trajectory.
trackStates()) {
507 if (pseudostate ==
nullptr) {
518 if ((pseudostate ==
nullptr) || pseudostate->fitQuality().chiSquared() < 10) {
523 std::unique_ptr<const TrackParameters> updpar(
m_updator->removeFromState(
525 pseudostate->measurement()->localParameters(),
526 pseudostate->measurement()->localCovariance()
529 if (updpar ==
nullptr) {
536 std::unique_ptr<const PseudoMeasurementOnTrack> newpseudo = std::make_unique<const PseudoMeasurementOnTrack>(
544 pseudostate->setMeasurement(std::move(newpseudo));
548 pseudostate->setMeasurementErrors(
errors);
549 pseudoupdated =
true;
560 *track->perigeeParameters(),
571 if (track !=
nullptr) {
572 track->info().addPatternReco(intrk1.
info());
573 track->info().addPatternReco(intrk2.
info());
584 const EventContext& ctx,
586 const Track & intrk1,
587 const Track & intrk2,
589 std::vector<MaterialEffectsOnTrack> & calomeots
591 ATH_MSG_DEBUG(
"--> entering GlobalChi2Fitter::mainCombinationStrategy");
596 const Track *indettrack = firstismuon ? &intrk2 : &intrk1;
597 const Track *muontrack = firstismuon ? &intrk1 : &intrk2;
599 auto [tmpfirstidpar, tmplastidpar] = getFirstLastIdPar(*indettrack);
600 std::unique_ptr<const TrackParameters> firstidpar =
unique_clone(tmpfirstidpar);
601 std::unique_ptr<const TrackParameters> lastidpar =
unique_clone(tmplastidpar);
603 if ((firstidpar ==
nullptr) || (lastidpar ==
nullptr)) {
617 std::unique_ptr<const TrackParameters> tp_closestmuon =
nullptr;
619 while (closestmuonmeas ==
nullptr) {
620 closestmuonmeas =
nullptr;
623 if ((**tsosit).measurementOnTrack() !=
nullptr) {
624 closestmuonmeas = (**tsosit).measurementOnTrack();
626 if (thispar !=
nullptr) {
627 const AmgVector(5) & parvec = thispar->parameters();
629 parvec[0], parvec[1], parvec[2], parvec[3], parvec[4], std::nullopt
643 std::unique_ptr<const TrackParameters> tmppar;
646 if ((tp_closestmuon !=
nullptr) && msEntranceIsValid) {
651 std::unique_ptr<const std::vector<const TrackStateOnSurface *>> matvec;
653 if (tmppar !=
nullptr) {
654 const Surface & associatedSurface = tmppar->associatedSurface();
655 std::unique_ptr<Surface> muonsurf =
nullptr;
661 double radius = cylbounds->
r();
663 muonsurf = std::make_unique<CylinderSurface>(trans,
radius + 1, hlength);
668 associatedSurface.
center().z() > 0 ?
669 associatedSurface.
center().z() + 1 :
670 associatedSurface.
center().z() - 1
674 associatedSurface.
center().x(),
675 associatedSurface.
center().y(),
679 trans.translation() << newpos;
682 double rmin = discbounds->
rMin();
683 double rmax = discbounds->
rMax();
684 muonsurf = std::make_unique<DiscSurface>(trans, rmin, rmax);
688 if (muonsurf !=
nullptr) {
700 std::vector<const TrackStateOnSurface *> tmp_matvec;
702 if ((matvec !=
nullptr) && !matvec->empty()) {
703 tmp_matvec = *matvec;
704 delete tmp_matvec.back();
705 tmp_matvec.pop_back();
707 for (
auto &
i : tmp_matvec) {
726 if (tmppar ==
nullptr) {
740 if (tmppar ==
nullptr) {
744 AmgVector(5) newpars = tmppar->parameters();
750 double newp2 = oldp * oldp + (!firstismuon ? 2 : -2) * de * std::sqrt(
mass *
mass + oldp * oldp) + de * de;
757 tp_closestmuon=tmppar->associatedSurface().createUniqueTrackParameters(
758 newpars[0], newpars[1], newpars[2], newpars[3], newpars[4], std::nullopt
774 for (; itStates != endState; ++itStates) {
781 if ((*itStates)->materialEffectsOnTrack() !=
nullptr) {
788 const auto *
const pBaseMEOT = (*itStates)->materialEffectsOnTrack();
792 const auto *
const pMEOT =
static_cast<const MaterialEffectsOnTrack *
>((*itStates)->materialEffectsOnTrack());
793 if ((pMEOT->scatteringAngles() ==
nullptr) or (pMEOT->energyLoss() ==
nullptr)) {
807 trajectory.
trackStates().back()->setTrackParameters(
nullptr);
810 std::unique_ptr<const TrackParameters> firstscatpar;
811 std::unique_ptr<const TrackParameters> lastscatpar;
814 double newqoverpid = 0;
817 double de = std::abs(calomeots[1].energyLoss()->deltaE());
818 double sigmade = std::abs(calomeots[1].energyLoss()->sigmaDeltaE());
820 double pbefore = std::abs(1 / firstidpar->parameters()[
Trk::qOverP]);
821 double pafter = std::abs(1 / tp_closestmuon->parameters()[
Trk::qOverP]);
822 double elosspull = (pbefore - pafter - de) / sigmade;
824 if (std::abs(elosspull) > 10) {
825 if (elosspull > 10) {
826 newqoverpid = 1 / (de + pafter + 10 * sigmade);
828 newqoverpid = 1 / (de + pafter - 10 * sigmade);
831 if (tp_closestmuon->parameters()[
Trk::qOverP] * newqoverpid < 0) {
835 const AmgVector(5) & newpar = firstidpar->parameters();
836 firstidpar=firstidpar->associatedSurface().createUniqueTrackParameters(
837 newpar[0], newpar[1], newpar[2], newpar[3], newqoverpid, std::nullopt
845 if (lastidpar ==
nullptr) {
851 *(firstismuon ? tp_closestmuon.get() : lastidpar.get()),
852 calomeots[0].associatedSurface(),
859 if (firstscatpar ==
nullptr) {
865 *(firstismuon ? firstidpar : tp_closestmuon),
866 calomeots[2].associatedSurface(),
873 if (lastscatpar ==
nullptr) {
877 std::optional<TransportJacobian> jac1, 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 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 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 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 double pull1 = std::abs(firstscatphi / firstscatmeff->sigmaDeltaPhi());
1109 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 std::unique_ptr<const TrackParameters> tmppar(
1311 calomeots[1].associatedSurface(),
1324 double oldp = std::abs(1 / tmppar->parameters()[
Trk::qOverP]);
1325 double de = std::abs(calomeots[1].energyLoss()->deltaE());
1327 double newp2 = oldp * oldp - 2 * de * std::sqrt(
mass *
mass + oldp * oldp) + de * de;
1333 double newqoverp =
sign / std::sqrt(newp2);
1338 tmppar->associatedSurface().createUniqueTrackParameters(
1345 calomeots[2].associatedSurface(),
1358 calomeots[2].associatedSurface(),
1372 calomeots[1].associatedSurface(),
1384 double newqoverp =
sign /
1385 (1. / std::abs(elosspar->parameters()[
Trk::qOverP]) +
1386 std::abs(calomeots[1].energyLoss()->deltaE()));
1390 std::unique_ptr<const TrackParameters>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) && (
1496 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 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());
1722 Track *track =
nullptr;
1730 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));
1787 track =
myfit(ctx, cache, trajectory, *firstidpar,
false,
muon);
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;
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)
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) {
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 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 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) {
2160 track.reset(
myfit(ctx, cache, trajectory, *track->perigeeParameters(),
false,
muon));
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) {
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>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 double sigmascat = std::sqrt(
m_scattool->sigmaSquare(
2612 newmeff = std::make_unique<GXFMaterialEffects>(newmeot);
2628 ((**trajectory.
trackStates().rbegin()).trackParameters() !=
nullptr)
2630 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);
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 double r = (std::abs(currentqoverp) > 1
e-10) ? -sintheta / (currentqoverp * 300. *
field[2]) : 1e6;
2973 double xc = parforextrap.
position().x() -
r * sinphi;
2974 double yc = parforextrap.
position().y() +
r * cosphi;
2977 double delta_s = (surf.
center().z() -
z0) / costheta;
2989 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 double r = (std::abs(currentqoverp) > 1
e-10) ? -sintheta / (currentqoverp * 300. *
field[2]) : 1e6;
3018 double xc = parforextrap.
position().x() -
r * sinphi;
3019 double yc = parforextrap.
position().y() +
r * cosphi;
3022 double d = xc * xc + yc * yc;
3023 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 double x1 = firstterm + secondterm;
3034 double x2 = firstterm - secondterm;
3035 firstterm = yc / 2 + (yc * (rcyl * rcyl -
r *
r)) / (2 *
d);
3036 secondterm = (mysqrt * xc) / (2 *
d);
3037 double y1 = firstterm - secondterm;
3038 double y2 = firstterm + secondterm;
3041 double dist1 = (
x -
x1) * (
x -
x1) + (
y -
y1) * (
y -
y1);
3042 double dist2 = (
x -
x2) * (
x -
x2) + (
y -
y2) * (
y -
y2);
3044 if (dist1 < dist2) {
3052 double phi1 = std::atan2(
y - yc,
x - xc);
3055 double delta_z =
r * deltaphi / tantheta;
3056 double z =
z0 + delta_z;
3069 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 double rmeas = oldstates[
i]->position().perp();
3107 double zmeas = oldstates[
i]->position().z();
3115 while (layerindex < (
int)
layers.size()) {
3117 double costracksurf = 0.0;
3138 if (oldstates[
i]->trackParameters() !=
nullptr) {
3139 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 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 double actualx0 =
X0 / costracksurf;
3205 double de = -std::abs(
3209 (
m_p != 0.0 ? std::abs(
m_p) : std::abs(1. / currentqoverp)),
3212 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 double lastz = laststate->
position().z();
3318 double lastr = laststate->
position().perp();
3325 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 = (
const DiscBounds *) (&(*it)->surfaceRepresentation().bounds());
3368 if (discbounds->
rMax() < firstr || discbounds->
rMin() > lastr) {
3372 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 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, lastmatindex = 0;
3522 std::vector<std::unique_ptr<GXFTrackState>> & oldstates = trajectory.
trackStates();
3536 for (
int i = 0;
i < (
int) oldstates.size();
i++) {
3537 if (oldstates[
i]->materialEffects() !=
nullptr) {
3546 if (firstsistate ==
nullptr) {
3547 if (oldstates[
i]->trackParameters() ==
nullptr) {
3548 std::unique_ptr<const TrackParameters> tmppar(
m_propagator->propagateParameters(
3551 oldstates[
i]->associatedSurface(),
3558 if (tmppar ==
nullptr)
return;
3560 oldstates[
i]->setTrackParameters(std::move(tmppar));
3562 firstsistate = oldstates[
i].get();
3564 lastsistate = oldstates[
i].get();
3572 if (lastsistate ==
nullptr) {
3573 throw std::logic_error(
"No track state");
3583 std::unique_ptr<const TrackParameters> tmppar(
m_propagator->propagateParameters(
3592 if (tmppar ==
nullptr)
return;
3604 indexoffset = lastmatindex;
3619 std::vector<std::pair<const Layer *, const Layer *>>
layers;
3620 std::vector<std::pair<const Layer *, const Layer *>> & upstreamlayers = trajectory.
upstreamMaterialLayers();
3634 const EventContext& ctx,
3640 if (refpar2 ==
nullptr) {
3650 std::vector<std::unique_ptr<GXFTrackState>> &
states = trajectory.
trackStates();
3651 std::vector<std::unique_ptr<GXFTrackState>> matstates;
3652 std::unique_ptr< const std::vector < const TrackStateOnSurface *>,
3653 void (*)(
const std::vector<const TrackStateOnSurface *> *) >
3654 matvec(
nullptr,&Trk::GlobalChi2Fitter::Cache::objVectorDeleter<TrackStateOnSurface>);
3655 bool matvec_used=
false;
3656 std::unique_ptr<TrackParameters> startmatpar1;
3657 std::unique_ptr<TrackParameters> startmatpar2;
3668 int npseudomuon1 = 0;
3669 int npseudomuon2 = 0;
3671 for (
auto & state :
states) {
3677 if (firstidhit ==
nullptr) {
3686 if (firsthit ==
nullptr) {
3687 firsthit = state->measurement();
3689 if (
tp ==
nullptr) {
3699 if (
tp ==
nullptr) {
3703 state->setTrackParameters(std::unique_ptr<const TrackParameters>(
tp));
3709 lasthit = state->measurement();
3715 if (firstidhit ==
nullptr) {
3716 firstidhit = state->measurement();
3719 if ((firstidpar ==
nullptr) && (
tp !=
nullptr)) {
3723 lastidhit = state->measurement();
3724 if (
tp !=
nullptr) {
3729 if (firstsiliconpar ==
nullptr) {
3730 firstsiliconpar =
tp;
3732 lastsiliconpar =
tp;
3744 if (firstmuonhit ==
nullptr) {
3745 firstmuonhit = state->measurement();
3746 if (
tp !=
nullptr) {
3750 lastmuonhit = state->measurement();
3751 if (
tp !=
nullptr) {
3757 if (meff->
deltaE() == 0) {
3758 if (firstcalopar ==
nullptr) {
3759 firstcalopar = state->trackParameters();
3761 lastcalopar = state->trackParameters();
3763 if (firstmatpar ==
nullptr) {
3764 firstmatpar = state->trackParameters();
3769 std::unique_ptr<TrackParameters> refpar;
3770 AmgVector(5) newpars = refpar2->parameters();
3777 newpars[0], newpars[1], newpars[2], newpars[3], newpars[4], std::nullopt
3780 if (firstmatpar !=
nullptr) {
3785 if ((startmatpar1 ==
nullptr) || ((firstidhit !=
nullptr) && (firstmuonhit !=
nullptr))) {
3791 const AmgVector(5) & newpars = startmatpar2->parameters();
3795 startmatpar2 = startmatpar2->associatedSurface().createUniqueTrackParameters(
3796 newpars[0], newpars[1], newpars[2], newpars[3],
3797 sign / std::sqrt(oldp * oldp + 2 * 100 *
MeV * std::sqrt(oldp * oldp +
mass *
mass) + 100 *
MeV * 100 *
MeV),
3802 AmgVector(5) newpars = startmatpar1->parameters();
3805 startmatpar1 = startmatpar1->associatedSurface().createUniqueTrackParameters(
3806 newpars[0], newpars[1], newpars[2], newpars[3], newpars[4], std::nullopt
3809 newpars = startmatpar2->parameters();
3812 startmatpar2 = startmatpar2->associatedSurface().createUniqueTrackParameters(
3813 newpars[0], newpars[1], newpars[2], newpars[3], newpars[4], std::nullopt
3821 refpar->momentum().unit()
3827 ATH_MSG_DEBUG(
"Obtaining upstream layers from Extrapolator");
3830 std::unique_ptr<const TrackParameters> tmppar;
3832 if (firstmuonhit !=
nullptr) {
3834 if (caloEntranceIsValid) {
3841 if (tmppar !=
nullptr) {
3842 destsurf = &tmppar->associatedSurface();
3847 if (matvec_used) cache.
m_matTempStore.push_back( std::move(matvec) );
3852 false, matEffects) );
3855 if (matvec && !matvec->empty()) {
3856 for (
int i = (
int)matvec->size() - 1;
i > -1;
i--) {
3861 std::unique_ptr<GXFMaterialEffects> meff = std::make_unique<GXFMaterialEffects>(*meot);
3862 const TrackParameters * newpars = (*matvec)[
i]->trackParameters() !=
nullptr ? (*matvec)[
i]->trackParameters()->
clone() :
nullptr;
3863 meff->setSigmaDeltaE(0);
3864 matstates.push_back(std::make_unique<GXFTrackState>(
3866 std::unique_ptr<const TrackParameters>(newpars)
3879 refpar->momentum().unit()
3885 ATH_MSG_DEBUG(
"Obtaining downstream ID layers from Extrapolator");
3887 std::unique_ptr<const TrackParameters> tmppar;
3888 std::unique_ptr<Surface> calosurf;
3889 if (firstmuonhit !=
nullptr) {
3891 if (caloEntranceIsValid) {
3899 if (tmppar !=
nullptr) {
3903 cylcalosurf =
static_cast<const CylinderSurface *
>(&tmppar->associatedSurface());
3908 disccalosurf =
static_cast<const DiscSurface *
>(&tmppar->associatedSurface());
3910 if (cylcalosurf !=
nullptr) {
3913 double radius = cylbounds.
r();
3915 calosurf = std::make_unique<CylinderSurface>(trans,
radius - 1, hlength);
3916 }
else if (disccalosurf !=
nullptr) {
3918 disccalosurf->
center().z() > 0 ?
3919 disccalosurf->
center().z() - 1 :
3920 disccalosurf->
center().z() + 1
3924 disccalosurf->
center().x(),
3925 disccalosurf->
center().y(),
3930 trans.translation() << newpos;
3933 double rmin = discbounds->
rMin();
3934 double rmax = discbounds->
rMax();
3935 calosurf = std::make_unique<DiscSurface>(trans, rmin, rmax);
3937 destsurf = calosurf.release();
3941 if (matvec_used) cache.
m_matTempStore.push_back( std::move(matvec) );
3943 ctx, *startmatpar2, *destsurf,
alongMomentum,
false, matEffects));
3944 matvec_used =
false;
3946 if (matvec && !matvec->empty()) {
3947 for (
const auto &
i : *matvec) {
3953 std::unique_ptr<GXFMaterialEffects> meff = std::make_unique<GXFMaterialEffects>(*meot);
3959 meff->setDeltaE(-5);
3962 meff->setScatteringSigmas(0, 0);
3965 meff->setSigmaDeltaE(50);
3968 const TrackParameters * newparams =
i->trackParameters() !=
nullptr ?
i->trackParameters()->clone() :
nullptr;
3970 matstates.push_back(std::make_unique<GXFTrackState>(
3972 std::unique_ptr<const TrackParameters>(newparams)
3984 if (cache.
m_calomat && (firstmuonhit !=
nullptr) && (firstidhit !=
nullptr)) {
3987 std::vector<MaterialEffectsOnTrack> calomeots =
m_calotool->extrapolationSurfacesAndEffects(
3996 if (calomeots.empty()) {
4001 std::unique_ptr<const TrackParameters> prevtrackpars =
unique_clone(lastidpar);
4002 if (lasthit == lastmuonhit) {
4003 for (
int i = 0;
i < (
int) calomeots.size();
i++) {
4006 std::unique_ptr<const TrackParameters> layerpar(
m_propagator->propagateParameters(
4009 calomeots[
i].associatedSurface(),
4016 if (layerpar ==
nullptr) {
4021 std::unique_ptr<GXFMaterialEffects> meff = std::make_unique<GXFMaterialEffects>(calomeots[
i]);
4024 lastcalopar = layerpar.get();
4028 double qoverp = layerpar->parameters()[
Trk::qOverP];
4029 double qoverpbrem = 0;
4033 (firstmuonpar !=
nullptr) &&
4034 std::abs(firstmuonpar->parameters()[
Trk::qOverP]) > 1.e-9
4036 qoverpbrem = firstmuonpar->parameters()[
Trk::qOverP];
4038 double sign = (qoverp > 0) ? 1 : -1;
4039 qoverpbrem =
sign / (1 / std::abs(qoverp) - std::abs(calomeots[
i].energyLoss()->deltaE()));
4042 const AmgVector(5) & newpar = layerpar->parameters();
4044 layerpar = layerpar->associatedSurface().createUniqueTrackParameters(
4045 newpar[0], newpar[1], newpar[2], newpar[3], qoverpbrem, std::nullopt
4047 meff->setdelta_p(1000 * (qoverpbrem - qoverp));
4050 matstates.push_back(std::make_unique<GXFTrackState>(
4052 std::unique_ptr<const TrackParameters>(layerpar !=
nullptr ? layerpar->clone() :
nullptr)
4054 prevtrackpars = std::move(layerpar);
4059 firsthit == firstmuonhit &&
4063 for (
int i = 0;
i < (
int) calomeots.size();
i++) {
4065 std::unique_ptr<const TrackParameters> layerpar(
m_propagator->propagateParameters(
4068 calomeots[
i].associatedSurface(),
4075 if (layerpar ==
nullptr) {
4080 std::unique_ptr<GXFMaterialEffects> meff = std::make_unique<GXFMaterialEffects>(calomeots[
i]);
4089 double qoverpbrem = layerpar->parameters()[
Trk::qOverP];
4094 (lastmuonpar !=
nullptr) &&
4095 std::abs(lastmuonpar->parameters()[
Trk::qOverP]) > 1.e-9
4099 double sign = (qoverpbrem > 0) ? 1 : -1;
4100 qoverp =
sign / (1 / std::abs(qoverpbrem) + std::abs(calomeots[
i].energyLoss()->deltaE()));
4103 meff->setdelta_p(1000 * (qoverpbrem - qoverp));
4104 const AmgVector(5) & newpar = layerpar->parameters();
4106 prevtrackpars = layerpar->associatedSurface().createUniqueTrackParameters(
4107 newpar[0], newpar[1], newpar[2], newpar[3], qoverp, std::nullopt
4111 matstates.insert(matstates.begin(), std::make_unique<GXFTrackState>(std::move(meff), std::move(layerpar)));
4116 if (lasthit == lastmuonhit && cache.
m_extmat) {
4117 std::unique_ptr<const Trk::TrackParameters> muonpar1;
4119 if (lastcalopar !=
nullptr) {
4121 if (msEntranceIsValid) {
4129 if (muonpar1 !=
nullptr) {
4137 rot.col(2) = trackdir;
4139 trans.linear().matrix() << rot;
4140 trans.translation() << muonpar1->
position() - .1 * trackdir;
4143 std::unique_ptr<const TrackParameters> curvlinpar(
m_extrapolator->extrapolateDirectly(
4151 if (curvlinpar !=
nullptr) {
4152 muonpar1 = std::move(curvlinpar);
4156 muonpar1 = std::unique_ptr<TrackParameters>(lastcalopar->
clone());
4160 muonpar1 = std::unique_ptr<TrackParameters>(refpar->clone());
4165 if (muonpar1 !=
nullptr) {
4175 0) and (firstmuonhit !=
nullptr)) {
4187 std::abs(distsol.
first()) < std::abs(distsol.
second()) ?
4194 if (firstmuonpar !=
nullptr) {
4195 AmgVector(5) newpars = firstmuonpar->parameters();
4202 newpars[0], newpars[1], newpars[2], newpars[3], newpars[4], std::nullopt
4205 std::unique_ptr<const TrackParameters> tmppar(
m_propagator->propagateParameters(
4215 if (tmppar !=
nullptr) {
4216 muonpar1 = std::move(tmppar);
4222 ATH_MSG_DEBUG(
"Obtaining downstream layers from Extrapolator");
4223 if (matvec_used) cache.
m_matTempStore.push_back( std::move(matvec) );
4226 states.back()->associatedSurface(),
4230 matvec_used =
false;
4237 if (matvec && !matvec->empty()) {
4238 for (
int j = 0; j < (
int) matvec->size(); j++) {
4244 std::unique_ptr<GXFMaterialEffects> meff = std::make_unique<GXFMaterialEffects>(*meot);
4249 std::abs(meff->deltaE()) > 25 &&
4250 std::abs((*matvec)[j]->trackParameters()->position().z()) < 13000
4255 const TrackParameters * newparams = (*matvec)[j]->trackParameters() !=
nullptr ? (*matvec)[j]->trackParameters()->
clone() :
nullptr;
4257 matstates.push_back(std::make_unique<GXFTrackState>(
4259 std::unique_ptr<const TrackParameters>(newparams)
4269 if (firsthit == firstmuonhit && cache.
m_extmat && (firstcalopar !=
nullptr)) {
4270 std::unique_ptr<const Trk::TrackParameters> muonpar1;
4273 if (msEntranceIsValid) {
4281 if (muonpar1 !=
nullptr) {
4289 rot.col(2) = trackdir;
4291 trans.linear().matrix() << rot;
4292 trans.translation() << muonpar1->
position() - .1 * trackdir;
4295 std::unique_ptr<const TrackParameters> curvlinpar(
m_extrapolator->extrapolateDirectly(
4303 if (curvlinpar !=
nullptr) {
4304 muonpar1 = std::move(curvlinpar);
4308 muonpar1 = std::unique_ptr<const TrackParameters>(firstcalopar->
clone());
4314 if (muonpar1 !=
nullptr) {
4325 ATH_MSG_DEBUG(
"Collecting upstream muon material from extrapolator");
4326 if (matvec_used) cache.
m_matTempStore.push_back( std::move(matvec) );
4329 states[0]->associatedSurface(),
4333 matvec_used =
false;
4335 if (matvec && !matvec->empty()) {
4336 ATH_MSG_DEBUG(
"Retrieved " << matvec->size() <<
" material states");
4338 for (
int j = 0; j < (
int) matvec->size(); j++) {
4341 if (meb !=
nullptr) {
4346 std::unique_ptr<GXFMaterialEffects> meff = std::make_unique<GXFMaterialEffects>(*meot);
4351 std::abs(meff->deltaE()) > 25 &&
4352 std::abs((*matvec)[j]->trackParameters()->position().z()) < 13000
4358 (*matvec)[j]->trackParameters() !=
nullptr
4359 ? (*matvec)[j]->trackParameters()->
clone()
4362 matstates.insert(matstates.begin(), std::make_unique<GXFTrackState>(
4364 std::unique_ptr<const TrackParameters>(tmpparams)
4377 std::vector<std::unique_ptr<GXFTrackState>> & newstates =
states;
4378 std::vector<std::unique_ptr<GXFTrackState>> oldstates = std::move(newstates);
4381 newstates.reserve(oldstates.size() + matstates.size());
4384 int firstlayerno = -1;
4394 bool addlayer =
true;
4396 while (addlayer && layerno < (
int) matstates.size()) {
4398 const TrackParameters *layerpar = matstates[layerno]->trackParameters();
4400 DistanceSolution distsol = oldstates[
i]->associatedSurface().straightLineDistanceEstimate(
4413 double trackimpact = std::abs(-refpar->position().x() * sinphi + refpar->position().y() * cosphi);
4415 if (trackimpact > cylinderradius - 5 *
mm) {
4421 if (
i == (
int) oldstates.size() - 1) {
4430 ATH_MSG_DEBUG(
"Material state with excessive scattering, skipping it");
4441 if (firstlayerno < 0) {
4442 firstlayerno = layerno;
4447 if ((layerpar !=
nullptr) && matEffects !=
pion && matEffects !=
muon) {
4449 const Layer *lay =
nullptr;
4451 if (tvol !=
nullptr) {
4467 if (matvec_used) cache.
m_matTempStore.push_back( std::move(matvec) );
4492 std::unique_ptr<const TrackParameters> per(
m_extrapolator->extrapolate(
4493 Gaudi::Hive::currentContext(),param, tmppersf,
oppositeMomentum,
false, matEffects));
4495 if (per ==
nullptr) {
4496 ATH_MSG_DEBUG(
"Cannot make Perigee with starting parameters");
4500 if(std::abs(per->position().z())>5000.) {
4501 ATH_MSG_WARNING(
"Pathological perigee well outside of tracking detector!! Returning nullptr");
4509 const EventContext& ctx,
4516 ATH_MSG_DEBUG(
"--> entering GlobalChi2Fitter::myfit_helper");
4544 if (trajectory.
nDOF() < 0) {
4554 ATH_MSG_WARNING(
"Attempt to apply material corrections with q/p=0, reject track");
4568 std::unique_ptr<const TrackParameters> per =
makePerigee(cache, param, matEffects);
4584 ATH_MSG_WARNING(
"Tracking Geometry Service and/or Material Updator Tool not configured");
4595 addMaterial(ctx, cache, trajectory, per !=
nullptr ? per.get() : ¶m, matEffects);
4598 ctx, cache, trajectory, per !=
nullptr ? per.get() : ¶m, matEffects);
4610 for (
const auto & state : trajectory.
trackStates()) {
4614 state->materialEffects()->deltaE() == 0
4616 scatstate2 = state.get();
4621 scatstate = state.get();
4628 if ((scatstate ==
nullptr) || (scatstate2 ==
nullptr)) {
4629 throw std::logic_error(
"Invalid scatterer");
4636 trajectory.
trackStates()[nstates / 2 - 1]->position() +
4642 std::unique_ptr<const TrackParameters> nearestpar;
4643 double mindist = 99999;
4644 std::vector < GXFTrackState * >mymatvec;
4647 if ((*it).trackParameters() ==
nullptr) {
4653 (*it).trackParameters()->position(),
4654 (*it).trackParameters()->momentum().unit())
4663 (((*it).measurement() !=
nullptr) && insideid) || (
4664 ((*it).materialEffects() !=
nullptr) &&
4666 (*it).materialEffects()->deltaE() == 0 ||
4667 ((*it).materialEffects()->sigmaDeltaPhi() == 0 &&
4669 (*it).materialEffects()->deltaPhi() != 0
4673 double dist = ((*it).trackParameters()->position() -
vertex).
perp();
4674 if (dist < mindist) {
4682 if (((*it).materialEffects() !=
nullptr) &&
distance > 0) {
4683 mymatvec.push_back(
it.get());
4687 if (nearestpar ==
nullptr) {
4691 for (
auto & state : mymatvec) {
4693 const Surface &matsurf = state->associatedSurface();
4695 nearestpar->position(), nearestpar->momentum().unit());
4703 std::unique_ptr<const TrackParameters> tmppar(
m_propagator->propagateParameters(
4713 if (tmppar ==
nullptr) {
4725 if (tmppar ==
nullptr) {
4735 AmgVector(5) newpars = tmppar->parameters();
4737 if (state->materialEffects()->sigmaDeltaE() > 0) {
4738 newpars[
Trk::qOverP] += .001 * state->materialEffects()->delta_p();
4741 double de = std::abs(state->materialEffects()->deltaE());
4743 double newp2 = oldp * oldp - 2 * de * std::sqrt(
mass *
mass + oldp * oldp) + de * de;
4749 nearestpar = tmppar->associatedSurface().createUniqueTrackParameters(
4750 newpars[0], newpars[1], newpars[2], newpars[3], newpars[4], std::nullopt
4754 std::unique_ptr<Trk::TrackParameters> tmpPars(
m_propagator->propagateParameters(
4765 if (tmpPars !=
nullptr) {
4766 per.reset(
static_cast < const Perigee *
>(tmpPars.release()));
4771 double oldp = 1. / std::abs(per->parameters()[
Trk::qOverP]);
4773 double newp = std::sqrt(oldp * oldp + 2 * toteloss * std::sqrt(oldp * oldp +
mass *
mass) + toteloss * toteloss);
4777 per = per->associatedSurface().createUniqueTrackParameters(
4782 if (per ==
nullptr) {
4799 }
else if (per ==
nullptr) {
4817 per = per->associatedSurface().createUniqueTrackParameters(
4824 if (per !=
nullptr) {
4834 Eigen::MatrixXd a_inv;
4835 a.resize(nfitpar, nfitpar);
4840 derivPool.setZero();
4842 for (std::unique_ptr<GXFTrackState> & state : trajectory.
trackStates()) {
4843 if (state->materialEffects() !=
nullptr) {
4846 state->setDerivatives(derivPool);
4849 bool doderiv =
true;
4881 double redchi2 = (trajectory.
nDOF() > 0) ? trajectory.
chi2() / trajectory.
nDOF() : 0;
4882 double prevredchi2 = (trajectory.
nDOF() > 0) ? trajectory.
prevchi2() / trajectory.
nDOF() : 0;
4891 (redchi2 < prevredchi2 &&
4892 (redchi2 > prevredchi2 - 1 || redchi2 < 2)) ||
4893 nsihits + ntrthits == nhits
4898 if (
it != 1 || nsihits != 0 || trajectory.
nDOF() <= 0 || trajectory.
chi2() / trajectory.
nDOF() <= 3) {
4913 if (ntrtprechits+ntrttubehits) {
4914 phf =
float(ntrtprechits)/
float(ntrtprechits+ntrttubehits);
4916 if (phf<m_minphfcut && it>=3) {
4917 if ((ntrtprechits+ntrttubehits)>=15) {
4922 <<
" | nTRTPrecHits = " << ntrtprechits
4923 <<
" | nTRTTubeHits = " << ntrttubehits
4924 <<
" | nOutliers = "
4944 if (trajectory.
prefit() == 0) {
4948 Eigen::LLT < Eigen::MatrixXd > lltOfW(
a);
4949 if (lltOfW.info() == Eigen::Success) {
4953 int ncols =
a.cols();
4954 Amg::MatrixX weightInvAMG = Amg::MatrixX::Identity(ncols, ncols);
4955 a_inv = lltOfW.solve(weightInvAMG);
4974 if (traj != &trajectory) {
4979 finaltrajectory = traj;
4987 if (nperpars == 5) {
4988 for (
int i = 0;
i <
a.cols();
i++) {
4989 a_inv(4,
i) *= .001;
4990 a_inv(
i, 4) *= .001;
4994 int scatterPos = nperpars + 2 * nscat;
4995 for (
int bremno = 0; bremno < nbrem; bremno++, scatterPos++) {
4996 for (
int i = 0;
i <
a.cols();
i++) {
4997 a_inv(scatterPos,
i) *= .001;
4998 a_inv(
i, scatterPos) *= .001;
5005 for (
int i = 0;
i < nperparams;
i++) {
5006 for (
int j = 0; j < nperparams; j++) {
5007 (errmat) (j,
i) = a_inv(j,
i);
5012 (errmat) (4, 4) = 1
e-20;
5016 std::unique_ptr<const TrackParameters> measper(
5018 perpars[0], perpars[1], perpars[2], perpars[3], perpars[4], std::move(errmat)
5028 std::unique_ptr<Track> track =
nullptr;
5030 if (finaltrajectory->
prefit() > 0) {
5031 if (finaltrajectory != &trajectory) {
5033 delete finaltrajectory;
5039 track =
makeTrack(ctx,cache, *finaltrajectory, matEffects);
5052 (track !=
nullptr) && (
5053 track->fitQuality()->numberDoF() != 0 &&
5054 track->fitQuality()->chiSquared() / track->fitQuality()->numberDoF() >
cut
5057 track.reset(
nullptr);
5061 if (track ==
nullptr) {
5065 if (finaltrajectory != &trajectory) {
5066 delete finaltrajectory;
5069 return track.release();
5073 const EventContext& ctx,
5074 const Cache & cache,
5078 int & bremno_maxbrempull,
5083 std::vector<std::unique_ptr<GXFTrackState>> &
states = trajectory.
trackStates();
5103 const int nmeas = (
int)
res.size();
5113 const int nDOF = trajectory.
nDOF();
5114 const bool doNewPseudoMeasurements = (
5118 std::abs((trajectory.
prevchi2() - trajectory.
chi2()) / nDOF) < 15 &&
5120 (nperpars == 0 || nidhits > 0)
5131 double maxbrempull = -0.2;
5145 for (
int hitno = 0; hitno < (
int)
states.size(); hitno++) {
5146 std::unique_ptr<GXFTrackState> & state =
states[hitno];
5160 doNewPseudoMeasurements &&
5162 !state->associatedSurface().isFree() &&
5163 !state->isRecalibrated()
5168 std::unique_ptr<const PseudoMeasurementOnTrack> newpseudo = std::make_unique<const PseudoMeasurementOnTrack>(
5174 state->setMeasurement(std::move(newpseudo));
5175 measbase = state->measurement();
5182 double *
errors = state->measurementErrors();
5184 for (
int i = 0;
i < 5;
i++) {
5204 res[measno] = residuals[
i];
5219 double *
errors = state->measurementErrors();
5220 for (
int i = 0;
i < 5;
i++) {
5233 ((trajectory.
prefit() == 0) || state->materialEffects()->deltaE() == 0)
5235 const double deltaPhi = state->materialEffects()->deltaPhi();
5236 const double measDeltaPhi = state->materialEffects()->measuredDeltaPhi();
5237 const double sigma2deltaPhi =
std::pow(state->materialEffects()->sigmaDeltaPhi(), 2);
5238 const double deltaTheta = state->materialEffects()->deltaTheta();
5239 const double sigma2deltaTheta =
std::pow(state->materialEffects()->sigmaDeltaTheta(), 2);
5241 if (trajectory.
prefit() != 1) {
5242 b[nperpars + 2 * scatno] -= (
deltaPhi - measDeltaPhi) / sigma2deltaPhi;
5243 b[nperpars + 2 * scatno + 1] -= deltaTheta / sigma2deltaTheta;
5245 b[nperpars + scatno] -= deltaTheta / sigma2deltaTheta;
5250 deltaTheta * deltaTheta / sigma2deltaTheta
5259 if ((state->materialEffects() !=
nullptr) && state->materialEffects()->sigmaDeltaE() > 0) {
5260 double averagenergyloss = std::abs(state->materialEffects()->deltaE());
5262 const double qoverp = limitInversePValue(qoverpbrem - state->materialEffects()->delta_p());
5263 const double pbrem = 1. / std::abs(qoverpbrem);
5264 const double p = 1. / std::abs(qoverp);
5265 const double mass = .001 * trajectory.
mass();
5267 const double bremEnergy = std::sqrt(pbrem * pbrem +
mass *
mass);
5269 const double resMaterial = .001 * averagenergyloss -
energy + bremEnergy;
5270 res[nmeas - nbrem + bremno] = resMaterial;
5272 const double sigde = state->materialEffects()->sigmaDeltaE();
5273 const double sigdepos = state->materialEffects()->sigmaDeltaEPos();
5274 const double sigdeneg = state->materialEffects()->sigmaDeltaENeg();
5276 double errorMaterial = .001 * state->materialEffects()->sigmaDeltaE();
5277 error[nmeas - nbrem + bremno] = errorMaterial;
5288 if (state->materialEffects()->isKink()) {
5289 maxbrempull = -999999999;
5290 state_maxbrempull =
nullptr;
5296 trajectory.
prefit() == 0 &&
5298 sigde != sigdepos &&
5301 const double elosspull = resMaterial / errorMaterial;
5303 if (trajectory.
mass() > 100) {
5308 if (std::abs(elosspull) > 1) {
5309 if (elosspull < -1) {
5310 state->materialEffects()->setSigmaDeltaE(sigdepos);
5312 state->materialEffects()->setSigmaDeltaE(sigdeneg);
5315 errorMaterial = .001 * state->materialEffects()->sigmaDeltaE();
5316 error[nmeas - nbrem + bremno] = errorMaterial;
5328 !state->materialEffects()->isKink() && (
5329 (
m_fixbrem == -1 && elosspull < maxbrempull) ||
5333 bremno_maxbrempull = bremno;
5334 state_maxbrempull = state.get();
5335 maxbrempull = elosspull;
5344 trajectory.
prefit() == 0 &&
5345 state->materialEffects()->sigmaDeltaPhi() == 0 &&
5346 state->materialEffects()->isMeasuredEloss() &&
5347 resMaterial / (.001 * state->materialEffects()->sigmaDeltaEAve()) > 2.5
5352 std::vector<MaterialEffectsOnTrack> calomeots =
5365 if (calomeots.size() == 3) {
5366 averagenergyloss = std::abs(calomeots[1].energyLoss()->deltaE());
5367 const double newres = .001 * averagenergyloss -
energy + bremEnergy;
5368 const double newerr = .001 * calomeots[1].energyLoss()->sigmaDeltaE();
5370 const double oldPull = resMaterial / errorMaterial;
5371 const double newPull = newres / newerr;
5373 if (std::abs(newPull) < std::abs(oldPull)) {
5374 ATH_MSG_DEBUG(
"Changing from measured to parametrized energy loss");
5376 state->materialEffects()->setEloss(std::unique_ptr<EnergyLoss>(calomeots[1].energyLoss()->
clone()));
5377 state->materialEffects()->setSigmaDeltaE(calomeots[1].energyLoss()->sigmaDeltaE());
5378 res[nmeas - nbrem + bremno] = newres;
5379 error[nmeas - nbrem + bremno] = newerr;
5383 state->materialEffects()->setMeasuredEloss(
false);
5393 for (
int imeas = 0; imeas < nmeas; imeas++) {
5394 if (
error[imeas] == 0) {
5409 const Cache & cache,
5415 const double oldChi2 = trajectory.
prevchi2();
5416 const double newChi2 = trajectory.
chi2();
5421 const double nDOF = trajectory.
nDOF();
5422 const double oldRedChi2 = (nDOF > 0) ? oldChi2 / nDOF : 0;
5423 const double newRedChi2 = (nDOF > 0) ? newChi2 / nDOF : 0;
5426 trajectory.
prefit() > 0 && (
5427 (newRedChi2 < 2 &&
it != 0) ||
5428 (newRedChi2 < oldRedChi2 + .1 && std::abs(newRedChi2 - oldRedChi2) < 1 &&
it != 1)
5441 int miniter = (nsihits != 0 && nsihits + ntrthits == nhits) ? 1 : 2;
5444 if (
it >= miniter && std::abs(oldChi2 - newChi2) < 1) {
5451 const int bremno_maxbrempull,
5457 if (state_maxbrempull ==
nullptr) {
5469 const int nmeas = (
int)
res.size();
5472 const double oldError =
error[nmeas - nbrem + bremno_maxbrempull];
5474 error[nmeas - nbrem + bremno_maxbrempull] = newError;
5477 if (
a.cols() != nFitPars) {
5481 const double errorRatio = oldError / newError;
5482 const double errorReductionRatio = 1 -
std::pow(errorRatio, 2);
5485 for (
int i = 0;
i < nFitPars;
i++) {
5486 if (weightderiv(nmeas - nbrem + bremno_maxbrempull,
i) == 0) {
5490 for (
int j =
i; j < nFitPars; j++) {
5491 const double newaij =
a(
i, j) - errorReductionRatio *
5492 weightderiv(nmeas - nbrem + bremno_maxbrempull,
i) *
5493 weightderiv(nmeas - nbrem + bremno_maxbrempull, j);
5495 a.fillSymmetric(
i, j, newaij);
5497 weightderiv(nmeas - nbrem + bremno_maxbrempull,
i) *= errorRatio;
5506 std::vector<std::unique_ptr<GXFTrackState>> &
states = trajectory.
trackStates();
5519 int nmeas = (
int) weightderiv.rows();
5523 for (std::unique_ptr<GXFTrackState> & state :
states) {
5527 const auto [scatmin, scatmax] = std::minmax(scatno, nscatupstream);
5528 const auto [bremmin, bremmax] = std::minmax(bremno, nbremupstream);
5535 const double sinStereo =
5537 state->sinStereo() :
5539 const double cosStereo =
5541 std::sqrt(1 -
std::pow(sinStereo, 2)) :
5549 auto getThisDeriv = [sinStereo, cosStereo, &derivatives](
int i,
int j) ->
double {
5550 if (
i == 0 && sinStereo != 0) {
5551 return derivatives(0, j) * cosStereo + sinStereo * derivatives(1, j);
5553 return derivatives(
i, j);
5557 for (
int i = 0;
i < 5;
i++) {
5572 if (
i == 0 && sinStereo != 0) {
5573 weightderiv.row(measno).head(
cols) =
5574 (derivatives.row(0).head(
cols) * cosStereo +
5575 sinStereo * derivatives.row(1).head(
cols)) /
5578 weightderiv.row(measno).head(
cols) = derivatives.row(
i).head(
cols) /
error[measno];
5582 for (
int j = scatmin; j < scatmax; j++) {
5583 if (trajectory.
prefit() == 1) {
5584 const int index = nperparams + j;
5587 const int index = nperparams + 2 * j;
5589 weightderiv(measno,
index + 1) = getThisDeriv(
i,
index + 1) /
error[measno];
5593 for (
int j = bremmin; j < bremmax; j++) {
5594 const int index = j + nperparams + 2 * nscat;
5601 double *
errors = state->measurementErrors();
5602 for (
int i = 0;
i < 5;
i++) {
5609 ((trajectory.
prefit() == 0) || state->materialEffects()->deltaE() == 0)
5614 if ((state->materialEffects() !=
nullptr) && state->materialEffects()->sigmaDeltaE() > 0) {
5616 double qoverpbrem = limitInversePValue(1000 * state->trackParameters()->parameters()[
Trk::qOverP]);
5617 double qoverp = limitInversePValue(qoverpbrem - state->materialEffects()->delta_p());
5619 double mass = .001 * trajectory.
mass();
5621 const auto thisMeasurementIdx{nmeas - nbrem + bremno};
5625 auto multiplier = [] (
double mass,
double qOverP){
5628 const auto qoverpTerm {multiplier(
mass, qoverp) /
error[thisMeasurementIdx]};
5629 const auto qoverpBremTerm {multiplier(
mass, qoverpbrem) /
error[thisMeasurementIdx]};
5632 weightderiv(thisMeasurementIdx, 4) = qoverpBremTerm - qoverpTerm;
5635 const auto bremNoBase = nperparams + 2 * nscat;
5636 if (bremno < nbremupstream) {
5637 weightderiv(thisMeasurementIdx, bremNoBase + bremno) = qoverpTerm;
5638 for (
int bremno2 = bremno + 1; bremno2 < nbremupstream; bremno2++) {
5639 weightderiv(thisMeasurementIdx, bremNoBase + bremno2) = qoverpTerm - qoverpBremTerm;
5642 weightderiv(thisMeasurementIdx, bremNoBase + bremno) = qoverpBremTerm;
5643 for (
int bremno2 = nbremupstream; bremno2 < bremno; bremno2++) {
5644 weightderiv(thisMeasurementIdx, bremNoBase + bremno2) = qoverpBremTerm - qoverpTerm;
5663 const int nMeas = (
int)
res.size();
5668 for (
int i = 0;
i < nPerPars;
i++) {
5677 const std::unique_ptr<GXFTrackState> & state = trajectory.
trackStates()[
i];
5680 if (meff ==
nullptr) {
5681 measno += state->numberOfMeasuredParameters();
5685 const int firstMeasurement =
i < nUpstreamStates ? 0 : measno;
5686 const int lastMeasurement =
i < nUpstreamStates ? measno : nMeas - nBrem;
5689 && (trajectory.
prefit() == 0 || meff->
deltaE() == 0)) {
5690 const int scatterPos = nPerPars + 2 * scatno;
5702 const int bremPos = nPerPars + nScatPars + bremno;
5713 const Cache & cache,
5726 const int nMeas = (
int)
res.size();
5728 for (
int k = 0;
k < nFitPars;
k++) {
5737 for (
int measno = minMeasK; measno < maxMeasK; measno++) {
5738 b[
k] +=
res[measno] * (1. /
error[measno]) * weightDeriv(measno,
k);
5748 if (
k == 4 ||
k >= nPerPars + nScatPars) {
5749 for (
int measno = nMeas - nBrem; measno < nMeas; measno++) {
5750 b[
k] +=
res[measno] * (1. /
error[measno]) * weightDeriv(measno,
k);
5757 const Cache & cache,
5764 for (
int k = 0;
k < nFitPars;
k++) {
5765 for (
int l =
k;
l < nFitPars;
l++) {
5770 for (
int measno = minMeas; measno < maxMeas; measno++) {
5771 a_kl += weightDeriv(measno,
k) * weightDeriv(measno,
l);
5774 a.fillSymmetric(
l,
k, a_kl);
5792 const int nMeas = (
int)
res.size();
5799 for (
int k = nPerPars;
k < nPerPars + nScatPars;
k += 2) {
5809 for (
int measno = nMeas - nBrem; measno < nMeas; measno++) {
5810 for (
int k = 4;
k < nFitPars;
k++) {
5812 k = nPerPars + nScatPars;
5815 for (
int l =
k;
l < nFitPars;
l++) {
5817 l = nPerPars + nScatPars;
5820 const double a_kl =
a(
l,
k) + weightDeriv(measno,
k) * weightDeriv(measno,
l);
5821 a.fillSymmetric(
l,
k, a_kl);
5833 const double oldRedChi2,
5834 const double newRedChi2
5842 bool weightChanged =
false;
5849 double newPhiWeight = 1.1;
5850 double newThetaWeight = 1.001;
5851 if (trajectory.
prefit() == 0) {
5857 newPhiWeight = 1.00000001;
5858 }
else if (
it == 1) {
5859 newPhiWeight = 1.0000001;
5860 }
else if (
it <= 3) {
5861 newPhiWeight = 1.0001;
5862 }
else if (
it <= 6) {
5863 newPhiWeight = 1.01;
5866 if (newRedChi2 > oldRedChi2 - 1 && newRedChi2 < oldRedChi2) {
5867 newPhiWeight = 1.0001;
5868 newThetaWeight = 1.0001;
5869 }
else if (newRedChi2 > oldRedChi2 - 25 && newRedChi2 < oldRedChi2) {
5870 newPhiWeight = 1.001;
5871 newThetaWeight = 1.0001;
5878 std::size_t scatno = 0;
5883 for (
const auto & state : trajectory.
trackStates()) {
5886 if (meff ==
nullptr) {
5890 const bool isValidPlaneSurface =
5892 static_cast<const PlaneSurface *
>(&state->associatedSurface()) !=
nullptr;
5897 if (meff->
deltaE() == 0 || (trajectory.
prefit() == 0 && isValidPlaneSurface)) {
5898 weightChanged =
true;
5900 const int scatNoIndex = 2 * scatno + nPerPars;
5905 message <<
"scatno is out of range " << scatno <<
" !< " << cache.
m_phiweight.size();
5906 throw std::range_error(
message.str());
5914 a(scatNoIndex, scatNoIndex) /= cache.
m_phiweight[scatno];
5918 a(scatNoIndex, scatNoIndex) *= newPhiWeight;
5919 }
else if (trajectory.
prefit() >= 2) {
5920 a(scatNoIndex, scatNoIndex) *= newPhiWeight;
5921 a(scatNoIndex + 1, scatNoIndex + 1) *= newThetaWeight;
5948 trajectory.
prefit() == 2 &&
5951 (newRedChi2 < oldRedChi2 - 25 || newRedChi2 > oldRedChi2)
5956 return weightChanged;
5965 std::size_t scatno = 0;
5976 message <<
"scatno is out of range " << scatno <<
" !< " << cache.
m_phiweight.size();
5977 throw std::range_error(
message.str());
5980 const bool isValidPlaneSurface =
5982 static_cast<const PlaneSurface *
>(&state->associatedSurface()) !=
nullptr;
5984 if (meff->
deltaE() == 0 || isValidPlaneSurface) {
5985 const int scatNoIndex = 2 * scatno + nPerPars;
5986 a(scatNoIndex, scatNoIndex) /= cache.
m_phiweight[scatno];
6001 const EventContext& ctx,
6010 const int nDOFold = trajectory.
nDOF();
6011 const double oldChi2 = trajectory.
chi2();
6012 const double oldRedChi2 = nDOFold > 0 ? oldChi2 / nDOFold : 0;
6034 int bremno_maxbrempull = 0;
6050 if ((state_maxbrempull !=
nullptr) && trajectory.
converged()) {
6059 const int nDOFnew = trajectory.
nDOF();
6060 const double newChi2 = trajectory.
chi2();
6061 const double newRedChi2 = nDOFnew > 0 ? newChi2 / nDOFnew : 0;
6063 ATH_MSG_DEBUG(
"old chi2: " << oldChi2 <<
"/" << nDOFold <<
"=" << oldRedChi2 <<
6064 ", new chi2: " << newChi2 <<
"/" << nDOFnew <<
"=" << newRedChi2);
6099 if (doDeriv || weightChanged) {
6110 if (trajectory.
prefit() == 0) {
6116 if (
nSiHits + nTrtHits != nHits) {
6123 (newRedChi2 < 2 || (newRedChi2 < oldRedChi2 && newRedChi2 > oldRedChi2 - .5))
6145 Eigen::LLT<Eigen::MatrixXd> llt(lu_m);
6147 if (llt.info() != Eigen::Success) {
6171 double d0 = refpar->parameters()[
Trk::d0];
6172 double z0 = refpar->parameters()[
Trk::z0];
6175 double qoverp = refpar->parameters()[
Trk::qOverP];
6177 if (nperparams > 0) {
6178 d0 += deltaParameters[0];
6179 z0 += deltaParameters[1];
6180 phi += deltaParameters[2];
6181 theta += deltaParameters[3];
6182 qoverp = (trajectory.
m_straightline) ? 0 : .001 * deltaParameters[4] + qoverp;
6194 std::vector < std::pair < double, double >>&scatangles = trajectory.
scatteringAngles();
6195 for (
int i = 0;
i < nscat;
i++) {
6196 scatangles[
i].first += deltaParameters[2 *
i + nperparams];
6197 scatangles[
i].second += deltaParameters[2 *
i + nperparams + 1];
6203 std::vector < double >&delta_ps = trajectory.
brems();
6204 for (
int i = 0;
i < nbrem;
i++) {
6205 delta_ps[
i] += deltaParameters[nperparams + 2 * nscat +
i];
6211 std::unique_ptr<const TrackParameters> newper(
6231 const EventContext& evtctx
6242 if (!splitProbContainer.
isValid()) {
6246 std::vector<std::unique_ptr<GXFTrackState>> &
states = trajectory.
trackStates();
6253 for (
size_t stateno = 0; stateno <
states.size(); stateno++) {
6258 measno +=
states[stateno-1]->numberOfMeasuredParameters();
6261 std::unique_ptr<GXFTrackState> & state =
states[stateno];
6273 const auto *
const rot =
static_cast<const RIO_OnTrack *
>(pMeas);
6285 if (!splitProb.isSplit()) {
6286 ATH_MSG_DEBUG(
"Pixel cluster is not split so no need to update" );
6290 std::unique_ptr < const RIO_OnTrack > newrot;
6291 double *olderror = state->measurementErrors();
6294 double newerror[5] = {-1,-1,-1,-1,-1};
6295 double newres[2] = {-1,-1};
6297 newrot.reset(
m_ROTcreator->correct(*prd, *trackpars, evtctx));
6304 newerror[0] = std::sqrt(
covmat(0, 0));
6305 newres[0] = newrot->localParameters()[
Trk::locX] - trackpars->parameters()[
Trk::locX];
6306 newerror[1] = std::sqrt(
covmat(1, 1));
6307 newres[1] = newrot->localParameters()[
Trk::locY] - trackpars->parameters()[
Trk::locY];
6309 if (
a.cols() != nfitpars) {
6314 for(
int k =0;
k<2;
k++ ){
6315 double oldres =
res[measno+
k];
6316 res[measno+
k] = newres[
k];
6317 err[measno+
k] = newerror[
k];
6319 for (
int i = 0;
i < nfitpars;
i++) {
6320 if (weightderiv(measno+
k,
i) == 0) {
6324 b[
i] -= weightderiv(measno+
k,
i) * (oldres / olderror[
k] - (newres[
k] * olderror[
k]) / (newerror[
k] * newerror[
k]));
6326 for (
int j =
i; j < nfitpars; j++) {
6330 weightderiv(measno+
k,
i) *
6331 weightderiv(measno+
k, j) *
6332 ((olderror[
k] * olderror[
k]) / (newerror[
k] * newerror[
k]) - 1)
6336 weightderiv(measno+
k,
i) *= olderror[
k] / newerror[
k];
6340 state->setMeasurement(std::move(newrot));
6341 state->setMeasurementErrors(newerror);
6356 const EventContext& ctx
6364 std::vector<std::unique_ptr<GXFTrackState>> &
states = trajectory.
trackStates();
6370 if (
a.cols() != nfitpars) {
6378 bool outlierremoved =
false;
6379 bool hitrecalibrated =
false;
6381 for (
int stateno = 0; stateno < (
int)
states.size(); stateno++) {
6382 std::unique_ptr<GXFTrackState> & state =
states[stateno];
6390 std::abs(state->trackParameters()->parameters()[
Trk::driftRadius]) > 1.05 * state->associatedSurface().bounds().r()
6395 outlierremoved =
true;
6397 double *
errors = state->measurementErrors();
6398 double olderror =
errors[0];
6402 for (
int i = 0;
i < nfitpars;
i++) {
6403 if (weightderiv(measno,
i) == 0) {
6407 b[
i] -=
res[measno] * weightderiv(measno,
i) / olderror;
6409 for (
int j =
i; j < nfitpars; j++) {
6412 a(
i, j) - weightderiv(measno,
i) * weightderiv(measno, j)
6415 weightderiv(measno,
i) = 0;
6419 }
else if (trtrecal) {
6420 double *
errors = state->measurementErrors();
6421 double olderror =
errors[0];
6423 const auto *
const thisMeasurement{state->measurement()};
6429 if (oldrot->prepRawData() !=
nullptr) {
6430 double dcradius = oldrot->prepRawData()->localPosition()[
Trk::driftRadius];
6432 double trackradius = state->trackParameters()->parameters()[
Trk::driftRadius];
6434 std::unique_ptr<const Trk::RIO_OnTrack> newrot =
nullptr;
6435 double distance = std::abs(std::abs(trackradius) - dcradius);
6437 if (
distance < scalefactor * dcerror && (olderror > 1. || trackradius * oldradius < 0)) {
6438 newrot.reset(
m_ROTcreator->correct(*oldrot->prepRawData(), *state->trackParameters(), ctx));
6439 }
else if (
distance > scalefactor * dcerror && olderror < 1.) {
6440 newrot.reset(
m_broadROTcreator->correct(*oldrot->prepRawData(), *state->trackParameters(), ctx));
6443 if (newrot !=
nullptr) {
6445 hitrecalibrated =
true;
6449 if ((measno < 0) or (measno >= (
int)
res.size())) {
6450 throw std::runtime_error(
6451 "'res' array index out of range in TrkGlobalChi2Fitter/src/GlobalChi2Fitter.cxx:" +
std::to_string(__LINE__)
6455 double oldres =
res[measno];
6456 double newres = newradius - state->trackParameters()->parameters()[
Trk::driftRadius];
6458 state->setMeasurement(std::move(newrot));
6462 for (
int i = 0;
i < nfitpars;
i++) {
6463 if (weightderiv(measno,
i) == 0) {
6467 b[
i] -= weightderiv(measno,
i) * (oldres / olderror - (newres * olderror) / (newerror * newerror));
6469 for (
int j =
i; j < nfitpars; j++) {
6476 i < nperpars + 2 * nscats &&
6477 (
i - nperpars) % 2 == 0
6484 a(
i, j) + weightderiv(measno,
i) * weightderiv(measno, j) * ((olderror * olderror) / (newerror * newerror) - 1) *
weight
6487 weightderiv(measno,
i) *= olderror / newerror;
6490 res[measno] = newres;
6491 err[measno] = newerror;
6500 measno += state->numberOfMeasuredParameters();
6504 if (trajectory.
nDOF() < 0) {
6509 if (outlierremoved || hitrecalibrated) {
6518 const EventContext& ctx,
6526 bool trackok =
false;
6528 std::unique_ptr < GXFTrajectory > cleanup_oldtrajectory;
6530 std::unique_ptr < GXFTrajectory > cleanup_newtrajectory;
6536 while (!trackok && oldtrajectory->
nDOF() > 0) {
6538 std::vector<std::unique_ptr<GXFTrackState>> &
states = oldtrajectory->
trackStates();
6546 if (nhits != nsihits) {
6550 double maxsipull = -1;
6552 int hitno_maxsipull = -1;
6553 int measno_maxsipull = -1;
6554 int stateno_maxsipull = 0;
6566 for (
int stateno = 0; stateno < (
int)
states.size(); stateno++) {
6567 std::unique_ptr<GXFTrackState> & state =
states[stateno];
6573 double *
errors = state->measurementErrors();
6575 const Amg::MatrixX & hitcov = state->measurement()->localCovariance();
6576 double sinstereo = state->sinStereo();
6577 double cosstereo = (sinstereo == 0) ? 1 : std::sqrt(1 - sinstereo * sinstereo);
6578 double weight1 = -1;
6580 if (hitcov(0, 0) > trackcov(0, 0)) {
6581 if (sinstereo == 0) {
6585 trackcov(0, 0) * cosstereo * cosstereo + 2 *
6586 trackcov(1, 0) * cosstereo * sinstereo + trackcov(1, 1) * sinstereo * sinstereo
6597 double sipull1 = weight1 > 0 ? std::abs(
res[measno] / std::sqrt(weight1)) : -1;
6600 std::abs(
res[measno + 1] / std::sqrt(weight2)) :
6603 sipull1 =
std::max(sipull1, sipull2);
6605 if (sipull1 > maxsipull) {
6606 maxsipull = sipull1;
6607 measno_maxsipull = measno;
6608 state_maxsipull = state.get();
6609 stateno_maxsipull = stateno;
6610 hitno_maxsipull = hitno;
6621 measno += state->numberOfMeasuredParameters();
6625 double maxpull = maxsipull;
6627 ATH_MSG_DEBUG(
" maxsipull: " << maxsipull <<
" hitno_maxsipull: " <<
6628 hitno_maxsipull <<
" n3sigma: " << n3sigma <<
" cut: " <<
cut <<
" cut2: " << cut2);
6639 state_maxsipull = oldtrajectory->
trackStates()[stateno_maxsipull].get();
6642 const auto *
const rot =
static_cast<const RIO_OnTrack *
>(pMeas);
6645 std::unique_ptr < const RIO_OnTrack > broadrot;
6650 Amg::VectorX parameterVector = trackpar_maxsipull->parameters();
6651 std::unique_ptr<const TrackParameters> trackparForCorrect(
6660 state_maxsipull->trackCovariance())
6664 newerror[0] = newerror[1] = newerror[2] = newerror[3] = newerror[4] = -1;
6665 double newpull = -1;
6666 double newpull1 = -1;
6667 double newpull2 = -1;
6668 double newres1 = -1;
6669 double newres2 = -1;
6670 double newsinstereo = 0;
6685 if (state_maxsipull->
sinStereo() != 0) {
6686 const auto [covEigenValueSmall, covStereoAngle] = principalComponentAnalysis2x2(
covmat);
6687 newerror[0] = std::sqrt(covEigenValueSmall);
6688 newsinstereo =
std::sin(covStereoAngle);
6690 newerror[0] = std::sqrt(
covmat(0, 0));
6693 double cosstereo = (newsinstereo == 0) ? 1. : std::sqrt(1 - newsinstereo * newsinstereo);
6695 if (cosstereo != 1.) {
6697 cosstereo * (broadrot->localParameters()[
Trk::locX] - trackpar_maxsipull->parameters()[
Trk::locX]) +
6698 newsinstereo * (broadrot->localParameters()[
Trk::locY] - trackpar_maxsipull->parameters()[
Trk::locY])
6701 newres1 = broadrot->localParameters()[
Trk::locX] - trackpar_maxsipull->parameters()[
Trk::locX];
6704 if (newerror[0] == 0.0) {
6705 ATH_MSG_WARNING(
"Measurement error is zero or negative, treating as outlier");
6708 newpull1 = std::abs(newres1 / newerror[0]);
6712 newerror[1] = std::sqrt(
covmat(1, 1));
6713 newres2 = broadrot->localParameters()[
Trk::locY] - trackpar_maxsipull->parameters()[
Trk::locY];
6714 newpull2 = std::abs(newres2 / newerror[1]);
6717 newpull =
std::max(newpull1, newpull2);
6723 (newerror[0] > 1.5 * olderror[0] || newerror[1] > 1.5 * std::abs(olderror[1]))
6725 if ((measno_maxsipull < 0) or(measno_maxsipull >= (
int)
res.size())) {
6726 throw std::runtime_error(
6727 "'res' array index out of range in TrkGlobalChi2Fitter/src/GlobalChi2Fitter.cxx:" +
std::to_string(__LINE__)
6732 newtrajectory = oldtrajectory;
6734 if (
a.cols() != nfitpars) {
6738 double oldres1 =
res[measno_maxsipull];
6739 res[measno_maxsipull] = newres1;
6740 err[measno_maxsipull] = newerror[0];
6742 for (
int i = 0;
i < nfitpars;
i++) {
6743 if (weightderiv(measno_maxsipull,
i) == 0) {
6747 b[
i] -= weightderiv(measno_maxsipull,
i) * (oldres1 / olderror[0] - (newres1 * olderror[0]) / (newerror[0] * newerror[0]));
6749 for (
int j =
i; j < nfitpars; j++) {
6753 weightderiv(measno_maxsipull,
i) *
6754 weightderiv(measno_maxsipull, j) *
6755 ((olderror[0] * olderror[0]) / (newerror[0] * newerror[0]) - 1)
6759 weightderiv(measno_maxsipull,
i) *= olderror[0] / newerror[0];
6763 double oldres2 =
res[measno_maxsipull + 1];
6764 res[measno_maxsipull + 1] = newres2;
6765 err[measno_maxsipull + 1] = newerror[1];
6767 for (
int i = 0;
i < nfitpars;
i++) {
6768 if (weightderiv(measno_maxsipull + 1,
i) == 0) {
6772 b[
i] -= weightderiv(measno_maxsipull + 1,
i) * (oldres2 / olderror[1] - (newres2 * olderror[1]) / (newerror[1] * newerror[1]));
6774 for (
int j =
i; j < nfitpars; j++) {
6778 weightderiv(measno_maxsipull + 1,
i) *
6779 weightderiv(measno_maxsipull + 1, j) *
6780 ((olderror[1] * olderror[1]) / (newerror[1] * newerror[1]) - 1)
6785 weightderiv(measno_maxsipull + 1,
i) *= olderror[1] / newerror[1];
6790 "Recovering outlier, hitno=" << hitno_maxsipull <<
" measno=" <<
6791 measno_maxsipull <<
" pull=" << maxsipull <<
" olderror_0=" <<
6792 olderror[0] <<
" newerror_0=" << newerror[0] <<
" olderror_1=" <<
6793 olderror[1] <<
" newerror_1=" << newerror[1]
6802 ((n3sigma < 2 && maxsipull > cut2 && maxsipull <
cut) || n3sigma > 1) &&
6812 "Removing outlier, hitno=" << hitno_maxsipull <<
", measno=" <<
6813 measno_maxsipull <<
" pull=" << maxsipull
6820 cleanup_newtrajectory = std::make_unique<GXFTrajectory>(*oldtrajectory);
6821 newtrajectory = cleanup_newtrajectory.get();
6823 if (newa.cols() != nfitpars) {
6829 if ((measno_maxsipull < 0) or(measno_maxsipull >= (
int)
res.size())) {
6830 throw std::runtime_error(
6831 "'res' array index out of range in TrkGlobalChi2Fitter/src/GlobalChi2Fitter.cxx:" +
std::to_string(__LINE__)
6835 double oldres1 =
res[measno_maxsipull];
6836 newres[measno_maxsipull] = 0;
6838 for (
int i = 0;
i < nfitpars;
i++) {
6839 if (weightderiv(measno_maxsipull,
i) == 0) {
6843 newb[
i] -= weightderiv(measno_maxsipull,
i) * oldres1 / olderror[0];
6845 for (
int j =
i; j < nfitpars; j++) {
6849 weightderiv(measno_maxsipull,
i) *
6850 weightderiv(measno_maxsipull, j)
6854 newweightderiv(measno_maxsipull,
i) = 0;
6858 double oldres2 =
res[measno_maxsipull + 1];
6859 newres[measno_maxsipull + 1] = 0;
6861 for (
int i = 0;
i < nfitpars;
i++) {
6862 if (weightderiv(measno_maxsipull + 1,
i) == 0) {
6866 newb[
i] -= weightderiv(measno_maxsipull + 1,
i) * oldres2 / olderror[1];
6868 for (
int j =
i; j < nfitpars; j++) {
6869 if (weightderiv(measno_maxsipull + 1, j) == 0) {
6876 weightderiv(measno_maxsipull + 1,
i) *
6877 weightderiv(measno_maxsipull + 1, j)
6881 newweightderiv(measno_maxsipull + 1,
i) = 0;
6885 newtrajectory->
setOutlier(stateno_maxsipull);
6909 ctx, cache, *newtrajectory,
it, *newap, *newbp, lu_m, doderiv);
6925 double oldchi2 = oldtrajectory->
chi2() / oldtrajectory->
nDOF();
6926 double newchi2 = (newtrajectory->
nDOF() > 0) ? newtrajectory->
chi2() / newtrajectory->
nDOF() : 0;
6929 if (newtrajectory->
nDOF() != oldtrajectory->
nDOF() && maxsipull > cut2) {
6930 mindiff = (oldchi2 > .33 *
m_chi2cut || noutl > 0) ? .8 : 1.;
6932 if (noutl == 0 && maxsipull <
cut - .5 && oldchi2 < .5 *
m_chi2cut) {
6937 if (newchi2 > oldchi2 || (newchi2 > oldchi2 - mindiff && newchi2 > .33 * oldchi2)) {
6938 ATH_MSG_DEBUG(
"Outlier not confirmed, keeping old trajectory");
6946 (void)cleanup_oldtrajectory.release();
6947 return oldtrajectory;
6949 if (oldtrajectory != newtrajectory) {
6950 cleanup_oldtrajectory = std::move(cleanup_newtrajectory);
6951 oldtrajectory = newtrajectory;
6958 Eigen::LLT < Eigen::MatrixXd > lltOfW(
a);
6959 if (lltOfW.info() == Eigen::Success) {
6963 int ncols =
a.cols();
6964 Amg::MatrixX weightInvAMG = Amg::MatrixX::Identity(ncols, ncols);
6965 fullcov = lltOfW.solve(weightInvAMG);
6983 oldtrajectory->
nDOF() > 0 &&
6992 (void)cleanup_oldtrajectory.release();
6993 return oldtrajectory;
7005 if (
const auto *pMeas{hit->measurement()};
7011 nrealmeas += hit->numberOfMeasuredParameters();
7021 if (
const auto *pMeas{hit->measurement()};
7027 for (
int i = measindex;
i < measindex + hit->numberOfMeasuredParameters();
i++) {
7030 if ((j == 4 && !oldtrajectory.
m_straightline) || j >= nperpars + 2 * nscat) {
7038 measindex += hit->numberOfMeasuredParameters();
7039 }
else if (hit->materialEffects() ==
nullptr) {
7040 measindex2 += hit->numberOfMeasuredParameters();
7046 const EventContext & ctx,
7051 GXFTrackState *firstmeasstate =
nullptr, *lastmeasstate =
nullptr;
7053 std::unique_ptr<const TrackParameters> per(
nullptr);
7056 std::unique_ptr<const TrackParameters> prevpar(
7061 std::vector<std::pair<const Layer *, const Layer *>> & upstreamlayers = oldtrajectory.
upstreamMaterialLayers();
7065 if (prevpar ==
nullptr) {
7070 const Layer *
layer = layer1 !=
nullptr ? layer1 : layer2;
7073 prevpar->position(), prevpar->momentum().unit()
7091 std::unique_ptr<const TrackParameters> layerpar(
7095 layer->surfaceRepresentation(),
7103 if (layerpar ==
nullptr) {
7107 if (
layer->surfaceRepresentation().bounds().inside(layerpar->localPosition())) {
7111 prevpar = std::move(layerpar);
7124 if (startfactor > 0.5) {
7125 std::unique_ptr<const TrackParameters> updatedpar =
m_matupdator->update(
7129 if (updatedpar !=
nullptr) {
7139 const Layer *endlayer = lastmeasstate->trackParameters()->associatedSurface().associatedLayer();
7149 if (endfactor > 0.5) {
7150 std::unique_ptr<const TrackParameters> updatedpar =
m_matupdator->update(
7151 lastmeasstate->trackParameters(), *endlayer,
alongMomentum, matEffects
7154 if (updatedpar !=
nullptr) {
7155 lastmeasstate->setTrackParameters(std::move(updatedpar));
7160 if (prevpar !=
nullptr) {
7172 if (per ==
nullptr) {
7173 ATH_MSG_DEBUG(
"Failed to extrapolate to perigee, returning 0");
7192 std::unique_ptr<GXFTrackState>
7194 const EventContext & ctx,
7201 if (per ==
nullptr) {
7205 ATH_MSG_DEBUG(
"Final perigee: " << *per <<
" pos: " << per->position() <<
" pT: " << per->pT());
7211 const std::vector<std::unique_ptr<TrackParameters>> & hc,
7212 std::set<Identifier> & id_set,
7213 std::set<Identifier> & sct_set,
7223 for (
const std::unique_ptr<TrackParameters> &
tp : hc) {
7229 if (
tp ==
nullptr) {
7240 if (de ==
nullptr) {
7251 if (id_set.find(
id) != id_set.end()) {
7300 if (sct_set.find(
os) != sct_set.end()) {
7301 ++
rv.m_sct_double_hole;
7330 for (
const std::unique_ptr<GXFTrackState> &
s : trajectory.
trackStates()) {
7342 std::vector<std::reference_wrapper<GXFTrackState>>
rv;
7350 for (
const std::unique_ptr<GXFTrackState> &
s : trajectory.
trackStates()) {
7364 rv.emplace_back(*
s);
7371 const TrkDetElementBase * de =
s->trackParameters()->associatedSurface().associatedDetectorElement();
7373 if (de !=
nullptr) {
7386 if (
s.get() == lastmeas) {
7396 const EventContext & ctx,
7397 const std::vector<std::reference_wrapper<GXFTrackState>> &
states
7409 constexpr
uint min_meas = 3;
7410 if (std::count_if(
states.begin(),
states.end(), [](
const GXFTrackState &
s){ return s.getStateType(TrackStateOnSurface::Measurement); }) < min_meas) {
7414 bool seen_meas =
false;
7416 std::set<Identifier> id_set;
7417 std::set<Identifier> sct_set;
7423 for (std::size_t
i = 0;
i <
states.size() - 1;
i++) {
7446 double dist = (
beg.trackParameters()->position() -
end.trackParameters()->position()).
norm();
7448 bool zStartValid = std::abs(
beg.trackParameters()->position().z())<10000.;
7450 ATH_MSG_DEBUG(
"Pathological track parameter well outside of detector");
7451 ATH_MSG_DEBUG(
"Propagator might have issue with this, skipping");
7461 if (seen_meas && dist >= 2.5 && zStartValid) {
7468 std::optional<std::vector<std::unique_ptr<TrackParameters>>> & hc =
beg.getHoles();
7469 std::vector<std::unique_ptr<TrackParameters>>
states;
7477 if (hc.has_value()) {
7516 std::vector<std::unique_ptr<Trk::TrackParameters>> bl =
m_extrapolator->extrapolateBlindly(
7537 const EventContext & ctx,
7543 auto trajectory = std::make_unique<Trk::TrackStates>();
7551 std::unique_ptr<GXFTrackState> perigee_ts =
makeTrackFindPerigee(ctx, cache, oldtrajectory, matEffects);
7553 if (perigee_ts ==
nullptr) {
7559 trajectory->reserve(tmptrajectory.
trackStates().size());
7565 hit->resetTrackCovariance();
7571 hit->materialEffects()))
7577 auto trackState = hit->trackStateOnSurface();
7578 hit->resetTrackCovariance();
7579 trajectory->emplace_back(trackState.release());
7582 auto qual = std::make_unique<FitQuality>(tmptrajectory.
chi2(), tmptrajectory.
nDOF());
7602 std::unique_ptr<Track>
rv = std::make_unique<Track>(
info, std::move(trajectory), std::move(
qual));
7612 std::unique_ptr<TrackSummary>
ts = std::make_unique<TrackSummary>();
7621 std::optional<TrackHoleCount> hole_count;
7646 if (hole_count.has_value()) {
7659 rv->setTrackSummary(std::move(
ts));
7665 GlobalChi2Fitter::~GlobalChi2Fitter() =
default;
7668 const EventContext & ctx,
7677 std::vector<std::unique_ptr<TrackParameters>>
rv =
m_extrapolator->extrapolateStepwise(
7691 &
rv.front()->associatedSurface() == &
src.associatedSurface() ||
7692 trackParametersClose(*
rv.front(),
src, 0.001) ||
7696 rv.front().reset(
nullptr);
7707 &
rv.back()->associatedSurface() == &
src.associatedSurface() ||
7708 trackParametersClose(*
rv.back(),
src, 0.001) ||
7712 rv.back().reset(
nullptr);
7719 const EventContext & ctx,
7727 std::unique_ptr<const TrackParameters>
rv;
7728 std::optional<TransportJacobian> jac{};
7739 if (
rv !=
nullptr && calcderiv) {
7744 std::optional<std::vector<std::unique_ptr<TrackParameters>>> extrapolation;
7753 std::move(extrapolation)
7758 const EventContext & ctx,
7769 ctx, prev,
ts, propdir, bf, calcderiv, holesearch
7772 if (
rv.m_parameters ==
nullptr) {
7773 propdir = invertPropdir(propdir);
7776 ctx, prev,
ts, propdir, bf, calcderiv, holesearch
7784 const EventContext& ctx,
7791 std::vector<std::unique_ptr<GXFTrackState>> &
states = trajectory.
trackStates();
7794 std::unique_ptr<const TrackParameters> tmptrackpar;
7796 for (
int hitno = nstatesupstream - 1; hitno >= 0; hitno--) {
7826 (
rv.m_parameters !=
nullptr) &&
7827 (prevtrackpar->
position() -
rv.m_parameters->position()).mag() > 5 *
mm
7833 if (
rv.m_parameters ==
nullptr) {
7834 ATH_MSG_DEBUG(
"propagation failed, prev par: " << *prevtrackpar <<
7835 " pos: " << prevtrackpar->
position() <<
" destination surface: " << surf1);
7839 states[hitno]->setTrackParameters(std::move(
rv.m_parameters));
7843 if (
rv.m_jacobian != std::nullopt) {
7845 states[hitno]->materialEffects() !=
nullptr &&
7846 states[hitno]->materialEffects()->deltaE() != 0 &&
7847 states[hitno]->materialEffects()->sigmaDeltaE() <= 0 &&
7850 double p = 1. / std::abs(currenttrackpar->parameters()[
Trk::qOverP]);
7851 double de = std::abs(
states[hitno]->materialEffects()->deltaE());
7853 double newp = std::sqrt(
p *
p + 2 * de * std::sqrt(
mass *
mass +
p *
p) + de * de);
7854 (*
rv.m_jacobian) (4, 4) = ((
p +
p * de / std::sqrt(
p *
p +
mass *
mass)) / newp) *
p *
p / (newp * newp);
7857 states[hitno]->setJacobian(*
rv.m_jacobian);
7858 }
else if (calcderiv) {
7865 if (meff !=
nullptr && hitno != 0) {
7867 surf, *meff, *
states[hitno]->trackParameters(), trajectory.
mass(), -1
7870 if (std::holds_alternative<FitterStatusCode>(
r)) {
7871 return std::get<FitterStatusCode>(
r);
7874 tmptrackpar = std::move(
std::get<std::unique_ptr<const TrackParameters>>(
r));
7875 prevtrackpar = tmptrackpar.get();
7877 prevtrackpar = currenttrackpar;
7883 for (
int hitno = nstatesupstream; hitno < (
int)
states.size(); hitno++) {
7905 (
rv.m_parameters !=
nullptr) &&
7907 (prevtrackpar->
position() -
rv.m_parameters->position()).mag() > 5 *
mm
7912 if (
rv.m_parameters ==
nullptr) {
7913 ATH_MSG_DEBUG(
"propagation failed, prev par: " << *prevtrackpar <<
7914 " pos: " << prevtrackpar->
7915 position() <<
" destination surface: " << surf);
7919 if (
rv.m_jacobian != std::nullopt) {
7921 states[hitno]->materialEffects() !=
nullptr &&
7922 states[hitno]->materialEffects()->deltaE() != 0 &&
7923 states[hitno]->materialEffects()->sigmaDeltaE() <= 0 &&
7926 double p = 1 / std::abs(
rv.m_parameters->parameters()[
Trk::qOverP]);
7927 double de = std::abs(
states[hitno]->materialEffects()->deltaE());
7929 double newp =
p *
p - 2 * de * std::sqrt(
mass *
mass +
p *
p) + de * de;
7932 newp = std::sqrt(newp);
7935 (*
rv.m_jacobian) (4, 4) = ((
p -
p * de / std::sqrt(
p *
p +
mass *
mass)) / newp) *
p *
p / (newp * newp);
7938 states[hitno]->setJacobian(*
rv.m_jacobian);
7939 }
else if (calcderiv) {
7946 if (meff !=
nullptr) {
7948 surf, *meff, *
rv.m_parameters, trajectory.
mass(), +1
7951 if (std::holds_alternative<FitterStatusCode>(
r)) {
7952 return std::get<FitterStatusCode>(
r);
7955 rv.m_parameters = std::move(
std::get<std::unique_ptr<const TrackParameters>>(
r));
7958 states[hitno]->setTrackParameters(std::move(
rv.m_parameters));
7959 prevtrackpar =
states[hitno]->trackParameters();
7978 ATH_MSG_DEBUG(
"Angles out of range, phi: " << newphi <<
" theta: " << newtheta);
7982 double newqoverp = 0;
7996 newqoverp = std::copysign(1 / std::sqrt(newp2),
old[
Trk::qOverP]);
8003 old[0],
old[1], newphi, newtheta, newqoverp, std::nullopt
8015 using Matrix55 = Eigen::Matrix<double, 5, 5>;
8017 Matrix55 initialjac;
8018 initialjac.setZero();
8019 initialjac(4, 4) = 1;
8021 Matrix55 jacvertex(initialjac);
8023 std::vector<Matrix55, Eigen::aligned_allocator<Matrix55>> jacscat(trajectory.
numberOfScatterers(), initialjac);
8024 std::vector<Matrix55, Eigen::aligned_allocator<Matrix55>> jacbrem(trajectory.
numberOfBrems(), initialjac);
8026 std::vector<std::unique_ptr<GXFTrackState>> &
states = trajectory.
trackStates();
8029 int hit_begin = 0, hit_end = 0, scatno = 0, bremno = 0;
8031 for (
bool forward : {
false,
true}) {
8033 hit_begin = nstatesupstream;
8035 scatno = nscatupstream;
8036 bremno = nbremupstream;
8038 hit_begin = nstatesupstream - 1;
8045 int hitno = hit_begin;
8046 forward ? (hitno < hit_end) : (hitno >= hit_end);
8047 hitno += (forward ? 1 : -1)
8050 state =
states[hitno].get();
8054 if (fillderivmat && state->derivatives().cols() != nfitpars) {
8055 state->derivatives().resize(5, nfitpars);
8056 state->derivatives().setZero();
8059 int jminscat = 0, jmaxscat = 4, jminbrem = 0, jmaxbrem = 4;
8061 if (hitno == (forward ? hit_end - 1 : 0)) {
8062 if (!fillderivmat) {
8070 Eigen::Matrix<double, 5, 5> & jac = state->jacobian();
8072 if (hitno == nstatesupstream + (forward ? 0 : -1)) {
8073 jacvertex.block<4, 5>(0, 0) = jac.block<4, 5>(0, 0);
8074 jacvertex(4, 4) = jac(4, 4);
8076 int jmin = 0, jmax = 0, jcnt = 0;
8077 int lp_bgn = 0, lp_end = 0;
8081 jcnt = jmax - jmin + 1;
8083 lp_bgn = forward ? nscatupstream : nscatupstream - 1;
8086 for (
int i = lp_bgn; forward ? (
i < lp_end) : (
i > lp_end);
i += (forward ? 1 : -1)) {
8088 i == scatno + (forward ? -1 : 1) &&
8089 prevstate !=
nullptr &&
8093 jacscat[
i].block(0, jmin, 4, jcnt) = jac.block(0, jmin, 4, jcnt);
8094 jacscat[
i](4, 4) = jac(4, 4);
8096 calculateJac(jac, jacscat[
i], jmin, jmax);
8100 Eigen::MatrixXd & derivmat = state->derivatives();
8101 int scatterPos = nperpars + 2 *
i;
8103 derivmat.block<4, 2>(0, scatterPos) = (forward ? 1 : -1) * jacscat[
i].block<4, 2>(0, 2);
8109 jcnt = jmax - jmin + 1;
8111 lp_bgn = forward ? nbremupstream : nbremupstream - 1;
8114 for (
int i = lp_bgn; forward ? (
i < lp_end) : (
i > lp_end);
i += (forward ? 1 : -1)) {
8116 i == bremno + (forward ? -1 : 1) &&
8121 jacbrem[
i].block(0, jmin, 4, jcnt) = jac.block(0, jmin, 4, jcnt);
8122 jacbrem[
i](4, 4) = jac(4, 4);
8124 calculateJac(jac, jacbrem[
i], jmin, jmax);
8128 Eigen::MatrixXd & derivmat = state->derivatives();
8129 int scatterPos = nperpars + 2 * nscats +
i;
8131 derivmat.block<5, 1>(0, scatterPos) = (forward ? .001 : -.001) * jacbrem[
i].block<5, 1>(0, 4);
8135 calculateJac(jac, jacvertex, 0, 4);
8139 Eigen::MatrixXd & derivmat = state->derivatives();
8140 derivmat.block(0, 0, 4, nperpars) = jacvertex.block(0, 0, 4, nperpars);
8142 if (nperpars == 5) {
8143 derivmat.col(4).segment(0, 4) *= .001;
8144 derivmat(4, 4) = .001 * jacvertex(4, 4);
8150 (!trajectory.
prefit() ||
states[hitno]->materialEffects()->deltaE() == 0)
8152 scatno += (forward ? 1 : -1);
8156 states[hitno]->materialEffects() &&
8157 states[hitno]->materialEffects()->sigmaDeltaE() > 0
8159 bremno += (forward ? 1 : -1);
8162 prevstate =
states[hitno].get();
8171 bool onlylocal)
const {
8177 std::vector<std::unique_ptr<GXFTrackState>> &
states = trajectory.
trackStates();
8181 int i = nstatesupstream;
8182 for (
int j = 0; j < (
int)
states.size(); j++) {
8183 if (j < nstatesupstream) {
8190 for (
int stateno = 0; stateno < (
int)
states.size(); stateno++) {
8191 if (stateno == 0 || stateno == nstatesupstream) {
8192 prevstate =
nullptr;
8195 std::unique_ptr<GXFTrackState> & state =
states[
index];
8196 if (state->materialEffects() !=
nullptr) {
8197 prevstate = state.get();
8201 if (!state->hasTrackCovariance()) {
8202 state->zeroTrackCovariance();
8204 AmgMatrix(5, 5) & trackerrmat = state->trackCovariance();
8206 if ((prevstate !=
nullptr) &&
8210 Eigen::Matrix<double, 5, 5> & jac = state->jacobian();
8213 trackerrmat = jac * prevcov * jac.transpose();
8217 trackerrmat = derivatives * fullcovmat * derivatives.transpose();
8228 bool errorok =
true;
8229 for (
int i = 0;
i < 5;
i++) {
8232 && trackerrmat(
i,
i) > meascov(j, j)) {
8234 double scale = std::sqrt(meascov(j, j) / trackerrmat(
i,
i));
8235 trackerrmat(
i,
i) = meascov(j, j);
8236 for (
int k = 0;
k < 5;
k++) {
8246 for (
int i = 0;
i < 5;
i++) {
8250 for (
int j = 0; j < 5; j++) {
8258 trackerrmat(4, 4) = 1
e-20;
8262 state->trackParameters();
8264 std::optional<
AmgMatrix(5, 5)> trkerrmat;
8266 if (state->hasTrackCovariance()) {
8267 trkerrmat = (state->trackCovariance());
8269 trkerrmat = std::nullopt;
8272 const AmgVector(5) & tpars = tmptrackpar->parameters();
8273 std::unique_ptr<const TrackParameters> trackpar(
8279 std::move(trkerrmat))
8281 state->setTrackParameters(std::move(trackpar));
8284 if (errorok && trajectory.
nDOF() > 0) {
8285 fitQual =
m_updator->fullStateFitQuality(
8286 *state->trackParameters(),
8294 state->setFitQuality(fitQual);
8296 prevstate = state.get();
8300 std::optional<TransportJacobian>
8302 const EventContext& ctx,
8316 std::optional<TransportJacobian> jac = std::make_optional<TransportJacobian>(J);
8319 0.01, 0.01, 0.00001, 0.00001, 0.000000001
8330 for (
int i = 0;
i < 5;
i++) {
8333 if (thisdiscsurf &&
i == 1) {
8337 vecpluseps[paraccessor.
pardef[
i]] += eps[
i];
8338 vecminuseps[paraccessor.
pardef[
i]] -= eps[
i];
8339 if (
i == 0 && thiscylsurf) {
8341 }
else if (
i == 1 && thisdiscsurf) {
8347 std::unique_ptr<const TrackParameters> parpluseps(
8357 std::unique_ptr<const TrackParameters> parminuseps(
8368 std::unique_ptr<const TrackParameters> newparpluseps(
8379 std::unique_ptr<const TrackParameters> newparminuseps(
8394 if (newparpluseps ==
nullptr) {
8406 if (newparminuseps ==
nullptr) {
8418 if ((newparpluseps ==
nullptr) || (newparminuseps ==
nullptr)) {
8422 for (
int j = 0; j < 5; j++) {
8423 double diff = newparpluseps->parameters()[paraccessor.
pardef[j]] -
8424 newparminuseps->parameters()[paraccessor.
pardef[j]];
8426 if (j == 0 && cylsurf) {
8428 }
else if (j == 1 && discsurf) {
8432 (*jac) (j,
i) =
diff / (2 * eps[
i]);
8445 (
"Configure the minimum number of Iterations via jobOptions");
8467 auto nmeas1 = pDataVector->size();
8468 const auto *pLastValue = (*pDataVector)[nmeas1 - 1];
8476 testrot =
static_cast<const RIO_OnTrack *
>(pLastValue);
8478 if (lastMeasIsCompetingRIO){
8484 if (testrot ==
nullptr) {
8485 const auto *pPenultimate = (*pDataVector)[nmeas1 - 2];
8488 if(penultimateIsRIO){
8489 testrot =
static_cast<const RIO_OnTrack *
>(pPenultimate);
8491 if (penultimateIsCompetingRIO){
8500 (testrot !=
nullptr) &&
8517 if (cond_obj ==
nullptr) {
8526 std::stringstream
msg;
8528 throw std::runtime_error(
msg.str());