56#include "CLHEP/Units/SystemOfUnits.h"
70#include <Eigen/StdVector>
81 return distsol.
first();
84 std::abs(distsol.
first()) < std::abs(distsol.
second()) ?
95 limitInversePValue(
double qOverP){
96 const double magnitude = std::abs(qOverP);
98 constexpr double maxP{100.*10e6*
MeV};
99 constexpr double minP{1.e-3*
MeV};
100 constexpr double lo {1./maxP};
101 constexpr double hi {1./minP};
102 const double limited = std::clamp(magnitude, lo, hi);
103 return std::copysign(limited, qOverP);
107 std::pair<const Trk::TrackParameters *, const Trk::TrackParameters *> getFirstLastIdPar(
const Trk::Track & track) {
113 while ((firstidpar ==
nullptr) && parit !=
track.trackParameters()->end()) {
115 ((**parit).covariance() !=
nullptr) &&
124 parit =
track.trackParameters()->end();
128 ((**parit).covariance() !=
nullptr) &&
133 }
while ((lastidpar ==
nullptr) && parit !=
track.trackParameters()->begin());
135 return std::make_pair(firstidpar, lastidpar);
150 std::abs(
a.parameters()[0] -
b.parameters()[0]) < e &&
151 std::abs(
a.parameters()[1] -
b.parameters()[1]) < e &&
152 std::abs(
a.parameters()[2] -
b.parameters()[2]) < e
163 calculateJac(Eigen::Matrix<double, 5, 5> &jac,
164 Eigen::Matrix<double, 5, 5> &out,
165 int jmin,
int jmax) {
169 out.block(0, 0, 4, jmin).setZero();
173 out.block(0, jmax + 1, 4, 5 - (jmax + 1)).setZero();
176 out(4, 4) = jac(4, 4);
186 std::pair<double, double> principalComponentAnalysis2x2(
const Amg::MatrixX & mat) {
187 const double trace =
mat(0, 0) +
mat(1, 1);
188 const double diagonalProduct =
mat(0, 0) *
mat(1, 1);
189 const double mat01Sq =
mat(0, 1) *
mat(0, 1);
190 const double discriminant = std::sqrt(trace * trace - 4. * (diagonalProduct - mat01Sq));
192 const double eigenValueSmall = 0.5 * (trace -
discriminant);
193 const double stereoAngle = 0.5 * std::asin(2 *
mat(0, 1) / (-discriminant));
195 return std::make_pair(eigenValueSmall, stereoAngle);
201 const std::string & t,
202 const std::string & n,
230 ATH_MSG_ERROR(
"Hole search requested but no boundary check tool provided.");
231 return StatusCode::FAILURE;
256 ATH_MSG_WARNING(
"FillDerivativeMatrix option selected, switching off acceleration!");
280 ATH_MSG_ERROR(
"Hole search requested but track summaries are disabled.");
281 return StatusCode::FAILURE;
286 return StatusCode::SUCCESS;
292 if (m_fit_status[
S_FITS] > 0) {
295 <<
" track fits failed because of a matrix inversion failure");
297 <<
" tracks were rejected by the outlier logic");
299 <<
" track fits failed because of a propagation failure");
301 <<
" track fits failed because of an invalid angle (theta/phi)");
303 <<
" track fits failed because the fit did not converge");
305 <<
" tracks did not pass the chi^2 cut");
307 <<
" tracks were killed by the energy loss update");
310 return StatusCode::SUCCESS;
315 std::unique_ptr<Track>
317 const EventContext& ctx,
323 ATH_MSG_DEBUG(
"--> entering GlobalChi2Fitter::fit(Track,Track,)");
338 const Track *indettrack = firstismuon ? &intrk2 : &intrk1;
339 const Track *muontrack = firstismuon ? &intrk1 : &intrk2;
341 bool measphi =
false;
348 rot = &crot->rioOnTrack(0);
358 const double dotprod2 = measdir.dot(
361 if (std::abs(dotprod1) < .5 && std::abs(dotprod2) < .5) {
369 auto [firstidpar, lastidpar] = getFirstLastIdPar(*indettrack);
371 if ((firstidpar ==
nullptr) || (lastidpar ==
nullptr)) {
375 std::unique_ptr<const TrackParameters> parforcalo =
unique_clone(firstismuon ? firstidpar : lastidpar);
378 const AmgVector(5) & newpars = parforcalo->parameters();
380 parforcalo=parforcalo->associatedSurface().createUniqueTrackParameters(
381 newpars[0], newpars[1], newpars[2], newpars[3], 1 / 5000., std::nullopt
385 std::vector < MaterialEffectsOnTrack > calomeots;
388 calomeots =
m_calotool->extrapolationSurfacesAndEffects(
392 parforcalo->associatedSurface(),
405 if (calomeots.empty()) {
410 std::unique_ptr<Track> track;
415 const bool tmp4 = cache.
m_idmat;
420 const double qoverpid = measperid !=
nullptr ? measperid->parameters()[
Trk::qOverP] : 0;
421 const double qoverpmuon = measpermuon !=
nullptr ? measpermuon->parameters()[
Trk::qOverP] : 0;
423 const AmgSymMatrix(5) * errmatid = measperid !=
nullptr ? measperid->covariance() :
nullptr;
424 const AmgSymMatrix(5) * errmatmuon = measpermuon !=
nullptr ? measpermuon->covariance() :
nullptr;
428 (errmatid !=
nullptr) &&
429 (errmatmuon !=
nullptr) &&
434 const double piderror = std::sqrt((*errmatid) (4, 4)) / (qoverpid * qoverpid);
435 const double pmuonerror = std::sqrt((*errmatmuon) (4, 4)) / (qoverpmuon * qoverpmuon);
436 const double energyerror = std::sqrt(
437 calomeots[1].energyLoss()->sigmaDeltaE() *
438 calomeots[1].energyLoss()->sigmaDeltaE() + piderror * piderror +
439 pmuonerror * pmuonerror
443 (std::abs(calomeots[1].energyLoss()->deltaE()) -
444 std::abs(1 / qoverpid) + std::abs(1 / qoverpmuon)
447 ATH_MSG_DEBUG(
"Changing from measured to parametrized energy loss");
452 parforcalo->associatedSurface(),
457 if (calomeots.empty()) {
465 bool firstfitwasattempted =
false;
468 if (!caloEntranceIsValid) {
479 qoverpid * qoverpmuon > 0
485 firstfitwasattempted =
true;
490 (track ==
nullptr) &&
491 !firstfitwasattempted &&
498 trajectory = trajectory2;
502 bool pseudoupdated =
false;
504 if (track !=
nullptr) {
505 for (std::unique_ptr<GXFTrackState> & pseudostate : trajectory.
trackStates()) {
506 if (pseudostate ==
nullptr) {
517 if ((pseudostate ==
nullptr) || pseudostate->fitQuality().chiSquared() < 10) {
522 const std::unique_ptr<const TrackParameters> updpar(
m_updator->removeFromState(
524 pseudostate->measurement()->localParameters(),
525 pseudostate->measurement()->localCovariance()
528 if (updpar ==
nullptr) {
533 covMatrix(0, 0) = 100;
535 std::unique_ptr<const PseudoMeasurementOnTrack> newpseudo = std::make_unique<const PseudoMeasurementOnTrack>(
539 std::move(covMatrix),
543 pseudostate->setMeasurement(std::move(newpseudo));
545 errors[0] = errors[2] = errors[3] = errors[4] = -1;
547 pseudostate->setMeasurementErrors(errors);
548 pseudoupdated =
true;
559 *track->perigeeParameters(),
570 if (track !=
nullptr) {
571 track->info().addPatternReco(intrk1.
info());
572 track->info().addPatternReco(intrk2.
info());
583 const EventContext& ctx,
585 const Track & intrk1,
586 const Track & intrk2,
588 std::vector<MaterialEffectsOnTrack> & calomeots
590 ATH_MSG_DEBUG(
"--> entering GlobalChi2Fitter::mainCombinationStrategy");
595 const Track *indettrack = firstismuon ? &intrk2 : &intrk1;
596 const Track *muontrack = firstismuon ? &intrk1 : &intrk2;
598 auto [tmpfirstidpar, tmplastidpar] = getFirstLastIdPar(*indettrack);
599 std::unique_ptr<const TrackParameters> firstidpar =
unique_clone(tmpfirstidpar);
600 std::unique_ptr<const TrackParameters> lastidpar =
unique_clone(tmplastidpar);
602 if ((firstidpar ==
nullptr) || (lastidpar ==
nullptr)) {
616 std::unique_ptr<const TrackParameters> tp_closestmuon =
nullptr;
618 while (closestmuonmeas ==
nullptr) {
619 closestmuonmeas =
nullptr;
622 if ((**tsosit).measurementOnTrack() !=
nullptr) {
623 closestmuonmeas = (**tsosit).measurementOnTrack();
625 if (thispar !=
nullptr) {
626 const AmgVector(5) & parvec = thispar->parameters();
628 parvec[0], parvec[1], parvec[2], parvec[3], parvec[4], std::nullopt
642 std::unique_ptr<const TrackParameters> tmppar;
645 if ((tp_closestmuon !=
nullptr) && msEntranceIsValid) {
650 std::unique_ptr<const std::vector<const TrackStateOnSurface *>> matvec;
652 if (tmppar !=
nullptr) {
653 const Surface & associatedSurface = tmppar->associatedSurface();
654 std::unique_ptr<Surface> muonsurf =
nullptr;
660 const double radius = cylbounds->
r();
662 muonsurf = std::make_unique<CylinderSurface>(trans, radius + 1, hlength);
666 const double newz = (
667 associatedSurface.
center().
z() > 0 ?
668 associatedSurface.
center().z() + 1 :
669 associatedSurface.
center().z() - 1
673 associatedSurface.
center().x(),
674 associatedSurface.
center().y(),
678 trans.translation() << newpos;
681 const double rmin = discbounds->
rMin();
682 const double rmax = discbounds->
rMax();
683 muonsurf = std::make_unique<DiscSurface>(trans, rmin, rmax);
687 if (muonsurf !=
nullptr) {
699 std::vector<const TrackStateOnSurface *> tmp_matvec;
701 if ((matvec !=
nullptr) && !matvec->empty()) {
702 tmp_matvec = *matvec;
703 delete tmp_matvec.back();
704 tmp_matvec.pop_back();
706 for (
auto & i : tmp_matvec) {
725 if (tmppar ==
nullptr) {
739 if (tmppar ==
nullptr) {
743 AmgVector(5) newpars = tmppar->parameters();
748 const double oldp = std::abs(1 / newpars[
Trk::qOverP]);
749 const double newp2 = oldp * oldp + (!firstismuon ? 2 : -2) * de * std::sqrt(mass * mass + oldp * oldp) + de * de;
756 tp_closestmuon=tmppar->associatedSurface().createUniqueTrackParameters(
757 newpars[0], newpars[1], newpars[2], newpars[3], newpars[4], std::nullopt
773 for (; itStates != endState; ++itStates) {
780 if ((*itStates)->materialEffectsOnTrack() !=
nullptr) {
787 const auto *
const pBaseMEOT = (*itStates)->materialEffectsOnTrack();
791 const auto *
const pMEOT =
static_cast<const MaterialEffectsOnTrack *
>((*itStates)->materialEffectsOnTrack());
792 if ((pMEOT->scatteringAngles() ==
nullptr) or (pMEOT->energyLoss() ==
nullptr)) {
806 trajectory.
trackStates().back()->setTrackParameters(
nullptr);
809 std::unique_ptr<const TrackParameters> firstscatpar;
810 std::unique_ptr<const TrackParameters> lastscatpar;
813 double newqoverpid = 0;
816 const double de = std::abs(calomeots[1].energyLoss()->deltaE());
817 const double sigmade = std::abs(calomeots[1].energyLoss()->sigmaDeltaE());
819 const double pbefore = std::abs(1 / firstidpar->parameters()[
Trk::qOverP]);
820 const double pafter = std::abs(1 / tp_closestmuon->parameters()[
Trk::qOverP]);
821 const double elosspull = (pbefore - pafter - de) / sigmade;
823 if (std::abs(elosspull) > 10) {
824 if (elosspull > 10) {
825 newqoverpid = 1 / (de + pafter + 10 * sigmade);
827 newqoverpid = 1 / (de + pafter - 10 * sigmade);
830 if (tp_closestmuon->parameters()[
Trk::qOverP] * newqoverpid < 0) {
834 const AmgVector(5) & newpar = firstidpar->parameters();
835 firstidpar=firstidpar->associatedSurface().createUniqueTrackParameters(
836 newpar[0], newpar[1], newpar[2], newpar[3], newqoverpid, std::nullopt
844 if (lastidpar ==
nullptr) {
850 *(firstismuon ? tp_closestmuon.get() : lastidpar.get()),
851 calomeots[0].associatedSurface(),
858 if (firstscatpar ==
nullptr) {
864 *(firstismuon ? firstidpar : tp_closestmuon),
865 calomeots[2].associatedSurface(),
872 if (lastscatpar ==
nullptr) {
876 std::optional<TransportJacobian> jac1;
877 std::optional<TransportJacobian> jac2;
878 std::unique_ptr<const TrackParameters> elosspar;
880 double firstscatphi = 0;
881 double secondscatphi = 0;
882 double firstscattheta = 0;
883 double secondscattheta = 0;
884 double muonscatphi = 0;
885 double muonscattheta = 0;
887 const TrackParameters *idscatpar = !firstismuon ? firstscatpar.get() : lastscatpar.get();
888 const TrackParameters *muonscatpar = firstismuon ? firstscatpar.get() : lastscatpar.get();
890 newqoverpid = idscatpar->parameters()[
Trk::qOverP];
892 const Amg::Vector3D calosegment = lastscatpar->position() - firstscatpar->position();
895 muonscattheta = calosegment.theta() - muonscatpar->parameters()[
Trk::theta];
898 for (
int i = 0; i < 2; i++) {
899 std::unique_ptr<const TrackParameters> tmpelosspar;
900 AmgVector(5) params1 = muonscatpar->parameters();
909 params1[0], params1[1], params1[2], params1[3], params1[4], std::nullopt
930 calomeots[1].associatedSurface(),
936 if ((tmpelosspar ==
nullptr) || (jac1 == std::nullopt)) {
945 const AmgVector(5) & newpars = tmpelosspar->parameters();
946 const std::unique_ptr<const TrackParameters> elosspar2(tmpelosspar->associatedSurface().createUniqueTrackParameters(
947 newpars[0], newpars[1], newpars[2], newpars[3], newqoverpid, std::nullopt
951 elosspar = std::move(tmpelosspar);
954 std::unique_ptr<const TrackParameters> scat2(
m_propagator->propagateParameters(
958 calomeots[0].associatedSurface() :
959 calomeots[2].associatedSurface(),
972 calomeots[0].associatedSurface() :
973 calomeots[2].associatedSurface(),
981 if ((scat2 ==
nullptr) || (jac2 == std::nullopt)) {
986 for (
int j = 0; j < 5; j++) {
987 for (
int k = 0; k < 5; k++) {
989 for (
int l = 0; l < 5; l++) {
990 jac3[j][k] += (*jac2) (j, l) * (*jac1) (l, k);
999 jac4(0, 0) = jac3[0][2];
1000 jac4(1, 1) = jac3[1][3];
1001 jac4(0, 1) = jac3[0][3];
1002 jac4(1, 0) = jac3[1][2];
1004 jac4 = jac4.inverse();
1016 discsurf =
static_cast<const Trk::DiscSurface *
>(&scat2->associatedSurface());
1018 if (cylsurf !=
nullptr) {
1019 dloc1 = -std::remainder(-dloc1, 2 *
M_PI * cylsurf->
bounds().
r());
1022 if (discsurf !=
nullptr) {
1023 dloc2 = -std::remainder(-dloc2, 2 *
M_PI);
1026 double dphi = jac4(0, 0) * dloc1 + jac4(0, 1) * dloc2;
1027 double dtheta = jac4(1, 0) * dloc1 + jac4(1, 1) * dloc2;
1033 muonscatphi += dphi;
1034 muonscattheta += dtheta;
1037 const double idscattheta = idscatpar->parameters()[
Trk::theta] - (scat2->parameters()[
Trk::theta] + dtheta);
1040 firstscatphi = muonscatphi;
1041 secondscatphi = idscatphi;
1042 firstscattheta = muonscattheta;
1043 secondscattheta = idscattheta;
1045 firstscatphi = -idscatphi;
1046 secondscatphi = -muonscatphi;
1047 firstscattheta = -idscattheta;
1048 secondscattheta = -muonscattheta;
1052 AmgVector(5) params2 = scat2->parameters();
1060 firstscatpar=scat2->associatedSurface().createUniqueTrackParameters(
1061 params2[0], params2[1], params2[2], params2[3], params2[4], std::nullopt
1063 idscatpar = firstscatpar.get();
1071 if (startPar !=
nullptr) {
1080 rot.col(2) = trackdir;
1083 trans.linear().matrix() << rot;
1084 trans.translation() << startPar->position() - .1 * trackdir;
1095 if (curvlinpar !=
nullptr) {
1096 startPar.reset(curvlinpar);
1100 firstscatpar = std::move(scat2);
1104 std::unique_ptr<GXFMaterialEffects> firstscatmeff = std::make_unique<GXFMaterialEffects>(calomeots[0]);
1105 std::unique_ptr<GXFMaterialEffects> elossmeff = std::make_unique<GXFMaterialEffects>(calomeots[1]);
1106 std::unique_ptr<GXFMaterialEffects> secondscatmeff = std::make_unique<GXFMaterialEffects>(calomeots[2]);
1108 const double pull1 = std::abs(firstscatphi / firstscatmeff->sigmaDeltaPhi());
1109 const double pull2 = std::abs(secondscatphi / secondscatmeff->sigmaDeltaPhi());
1112 for (
auto & i : tmp_matvec) {
1117 firstscatmeff->setScatteringAngles(firstscatphi, firstscattheta);
1118 secondscatmeff->setScatteringAngles(secondscatphi, secondscattheta);
1121 elossmeff->setdelta_p(1000 * (lastscatpar->parameters()[
Trk::qOverP] - newqoverpid));
1123 elossmeff->setdelta_p(1000 * (newqoverpid - firstscatpar->parameters()[
Trk::qOverP]));
1126 elossmeff->setSigmaDeltaE(calomeots[1].energyLoss()->sigmaDeltaE());
1128 trajectory.
addMaterialState(std::make_unique<GXFTrackState>(std::move(firstscatmeff), std::move(firstscatpar)), -1);
1129 trajectory.
addMaterialState(std::make_unique<GXFTrackState>(std::move(elossmeff), std::move(elosspar)), -1);
1130 trajectory.
addMaterialState(std::make_unique<GXFTrackState>(std::move(secondscatmeff), std::move(lastscatpar)), -1);
1133 for (
auto & i : tmp_matvec) {
1140 if (startPar ==
nullptr) {
1145 (pull1 > 5 || pull2 > 5) &&
1151 bool largegap =
false;
1152 double previousz = 0;
1154 for (itStates2 = beginStates2; itStates2 != endState2; ++itStates2) {
1159 if (meff !=
nullptr) {
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());
1196 covMatrix(0, 0) = 100;
1198 std::unique_ptr<const PseudoMeasurementOnTrack> newpseudo = std::make_unique<const PseudoMeasurementOnTrack>(
1200 std::move(covMatrix),
1201 pseudopar->associatedSurface()
1204 std::unique_ptr<GXFTrackState> pseudostate = std::make_unique<GXFTrackState>(std::move(newpseudo), std::move(pseudopar));
1208 errors[0] = errors[2] = errors[3] = errors[4] = -1;
1211 pseudostate->setMeasurementErrors(errors);
1218 std::abs(trajectory.
trackStates().back()->position().z()) > 20000 &&
1219 std::abs(previousz) < 12000
1225 previousz = trajectory.
trackStates().back()->position().z();
1243 const EventContext& ctx,
1245 const Track & intrk1,
1246 const Track & intrk2,
1248 std::vector<MaterialEffectsOnTrack> & calomeots
1250 ATH_MSG_DEBUG(
"--> entering GlobalChi2Fitter::backupCombinationStrategy");
1253 const Track *indettrack = firstismuon ? &intrk2 : &intrk1;
1266 if (pParametersVector->size() > 1)
1267 firstidpar = (*pParametersVector)[1];
1269 firstidpar = pParametersVector->back();
1271 std::unique_ptr<const TrackParameters> lastidpar =
nullptr;
1272 if ((firstidpar !=
nullptr) && (cache.
m_caloEntrance !=
nullptr))
1276 if (lastidpar ==
nullptr) {
1277 lastidpar.reset(pParametersVector->back()->clone());
1280 std::unique_ptr < const TrackParameters > firstscatpar;
1281 std::unique_ptr < const TrackParameters > lastscatpar;
1282 std::unique_ptr < const TrackParameters > elosspar;
1287 lastidpar->position(),
1288 lastidpar->momentum(),
1290 lastidpar->position()
1297 calomeots[0].associatedSurface(),
1303 if (!firstscatpar) {
1307 const std::unique_ptr<const TrackParameters> tmppar(
1311 calomeots[1].associatedSurface(),
1324 const double oldp = std::abs(1 / tmppar->parameters()[
Trk::qOverP]);
1325 const double de = std::abs(calomeots[1].energyLoss()->deltaE());
1327 double newp2 = oldp * oldp - 2 * de * std::sqrt(mass * mass + oldp * oldp) + de * de;
1333 const double newqoverp =
sign / std::sqrt(newp2);
1335 const AmgVector(5) & pars = tmppar->parameters();
1338 tmppar->associatedSurface().createUniqueTrackParameters(
1339 pars[0], pars[1], pars[2], pars[3], newqoverp, std::nullopt
1345 calomeots[2].associatedSurface(),
1358 calomeots[2].associatedSurface(),
1372 calomeots[1].associatedSurface(),
1383 const double sign = (elosspar->parameters()[
Trk::qOverP] < 0) ? -1 : 1;
1384 const double newqoverp =
sign /
1385 (1. / std::abs(elosspar->parameters()[
Trk::qOverP]) +
1386 std::abs(calomeots[1].energyLoss()->deltaE()));
1388 const AmgVector(5) & pars = elosspar->parameters();
1390 std::unique_ptr<const TrackParameters>
const tmppar(
1391 elosspar->associatedSurface().createUniqueTrackParameters(
1392 pars[0], pars[1], pars[2], pars[3], newqoverp, std::nullopt
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);
1482 if ((rot !=
nullptr) &&
m_DetID->is_mdt(rot->
identify()) && (triggersurf1 !=
nullptr)) {
1486 (rot !=
nullptr) && (
1493 const double dotprod1 = measdir.dot(
Amg::Vector3D(0, 0, 1));
1496 const bool measphi = std::abs(dotprod1) <= .5 && std::abs(dotprod2) <= .5;
1500 (*itStates2)->trackParameters() !=
nullptr ?
1501 (*itStates2)->trackParameters()->position() :
1503 if (triggersurf1 !=
nullptr) {
1504 triggerpos2 = thispos;
1505 triggersurf2 = surf;
1507 mdtbetweenphihits =
true;
1510 triggerpos1 = thispos;
1511 triggersurf1 = surf;
1517 double mdttrig1 = 999999;
1518 double mdttrig2 = 999999;
1519 const Surface *mdtsurf1 =
nullptr;
1520 const Surface *mdtsurf2 =
nullptr;
1523 itStates2 = (!firstismuon ? beginStates2 : endState - 1);
1524 itStates2 != (!firstismuon ? endState2 : beginStates - 1);
1525 (!firstismuon ? ++itStates2 : --itStates2)
1527 const Surface *surf =
nullptr;
1529 ((*itStates2)->measurementOnTrack() !=
nullptr) &&
1532 surf = &(*itStates2)->measurementOnTrack()->associatedSurface();
1535 if (surf ==
nullptr) {
1538 const auto *
const pThisMeasurement = (*itStates2)->measurementOnTrack();
1542 if (isCompetingRioOnTrack) {
1544 rot = &crot->rioOnTrack(0);
1547 rot =
static_cast<const RIO_OnTrack *
>(pThisMeasurement);
1553 (*itStates2)->trackParameters() !=
nullptr ?
1554 (*itStates2)->trackParameters()->position() :
1555 pThisMeasurement->globalPosition();
1556 if (triggerpos1.mag() > 1 && (globpos - triggerpos1).mag() < mdttrig1) {
1557 mdttrig1 = (globpos - triggerpos1).
mag();
1560 if (triggerpos2.mag() > 1 && (globpos - triggerpos2).mag() < mdttrig2) {
1561 mdttrig2 = (globpos - triggerpos2).
mag();
1568 std::vector<GXFTrackState *> outlierstates;
1569 std::vector<GXFTrackState *> outlierstates2;
1571 outlierstates.reserve(10);
1572 outlierstates2.reserve(10);
1574 std::unique_ptr<PseudoMeasurementOnTrack> newpseudo;
1576 for (itStates2 = beginStates2; itStates2 != endState2; ++itStates2) {
1577 const auto *
const pMeasurement{(*itStates2)->measurementOnTrack()};
1579 const bool isStraightLine =
1580 pMeasurement !=
nullptr ?
1587 (newpseudo ==
nullptr) && (
1588 (itStates2 == beginStates2 || itStates2 == beginStates2 + 1) &&
1589 std::abs(pMeasurement->globalPosition().z()) < 10000
1592 std::unique_ptr<const TrackParameters> par2;
1593 if (((*itStates2)->trackParameters() !=
nullptr) && nphi > 99) {
1594 par2.reset((*itStates2)->trackParameters()->clone());
1606 if (par2 ==
nullptr) {
1610 covMatrix(0, 0) = 100;
1612 newpseudo = std::make_unique<PseudoMeasurementOnTrack>(
1614 std::move(covMatrix),
1615 par2->associatedSurface()
1618 std::unique_ptr<GXFTrackState> firstpseudo = std::make_unique<GXFTrackState>(std::move(newpseudo), std::move(par2));
1622 errors[0] = errors[2] = errors[3] = errors[4] = -1;
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;
1655 covMatrix(0, 0) = 100;
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);
1665 errors[0] = errors[2] = errors[3] = errors[4] = -1;
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;
1683 covMatrix(0, 0) = 100;
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);
1693 errors[0] = errors[2] = errors[3] = errors[4] = -1;
1696 pseudostate2->setMeasurementErrors(errors);
1698 outlierstates2.push_back(pseudostate2.get());
1714 outlierstates.push_back(trajectory.
trackStates().back().get());
1722 Track *track =
nullptr;
1730 const std::unique_ptr<Trk::Track> tmp_track(
1731 myfit(ctx, cache, trajectory, *startpar2,
false,
muon));
1738 std::abs(trajectory.
residuals().tail<1>()(0) / trajectory.
errors().tail<1>()(0)) > 10
1744 if (firstpseudostate !=
nullptr) {
1747 covMatrix(0, 0) = 100;
1749 std::unique_ptr<const PseudoMeasurementOnTrack> newpseudo = std::make_unique<const PseudoMeasurementOnTrack>(
1751 std::move(covMatrix),
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 alignCache.m_derivMatrix.reset();
1830 alignCache.m_fullCovarianceMatrix.reset();
1831 alignCache.m_iterationsOfLastFit = 0;
1834 fitIm(ctx, cache, inputTrack, runOutlier, matEffects);
1835 if(newTrack !=
nullptr){
1837 alignCache.m_derivMatrix = std::make_unique<Amg::MatrixX>(cache.
m_derivmat);
1839 alignCache.m_fullCovarianceMatrix = std::make_unique<Amg::MatrixX>(cache.
m_fullcovmat);
1840 alignCache.m_iterationsOfLastFit = cache.
m_lastiter;
1847 const EventContext& ctx,
1849 const Track& inputTrack,
1854 ATH_MSG_DEBUG(
"--> entering GlobalChi2Fitter::fit(Track,,)");
1870 ATH_MSG_WARNING(
"Track without track parameters, cannot perform fit");
1878 if (minpar ==
nullptr) {
1896 const Surface *firsthitsurf =
nullptr;
1897 const Surface *lasthitsurf =
nullptr;
1899 bool hasmuon =
false;
1900 bool iscombined =
false;
1901 bool seenphimeas =
false;
1905 for (; itStates != endState; ++itStates) {
1906 const auto *
const pMeasurement = (**itStates).measurementOnTrack();
1908 (pMeasurement ==
nullptr) &&
1909 ((**itStates).materialEffectsOnTrack() !=
nullptr) &&
1915 if (pMeasurement !=
nullptr) {
1916 const Surface *surf = &pMeasurement->associatedSurface();
1925 if (firsthitsurf ==
nullptr) {
1926 firsthitsurf = surf;
1929 if (
m_DetID->is_indet(hitid)) {
1934 if ((**itStates).trackParameters() !=
nullptr) {
1935 lastidpar = (**itStates).trackParameters();
1936 if (firstidpar ==
nullptr) {
1937 firstidpar = lastidpar;
1943 const double dotprod1 = measdir.dot(
Amg::Vector3D(0, 0, 1));
1945 if (std::abs(dotprod1) < .5 && std::abs(dotprod2) < .5) {
1947 if (std::abs(surf->
center().z()) > 13000) {
1950 if (surf->
center().perp() > 9000 && std::abs(surf->
center().z()) < 13000) {
1962 if (iscombined && seenphimeas && (phiem || phibo)) {
1978 if (firstidpar == lastidpar) {
1979 firstidpar = lastidpar =
nullptr;
1989 (firstidpar !=
nullptr)
2000 const bool tmpsirecal = cache.
m_sirecal;
2001 std::unique_ptr<Track> tmptrack =
nullptr;
2012 tmptrack.reset(
myfit(ctx, cache, trajectory, *minpar,
false, matEffects));
2015 if (tmptrack ==
nullptr) {
2024 bool isbrem =
false;
2026 unsigned int n_brem=0;
2028 for (std::unique_ptr<GXFTrackState> & state : trajectory.
trackStates()) {
2031 if (meff !=
nullptr) {
2037 const double p = 1. / std::abs(layerpars->parameters()[
Trk::qOverP] - .0005 * meff->
delta_p());
2039 std::optional<Amg::Vector2D> locpos(state->associatedSurface().globalToLocal(layerpars->
position()));
2040 const Amg::Vector3D layerNormal(state->associatedSurface().normal(*locpos));
2041 double costracksurf = 1.;
2043 costracksurf = std::abs(layerNormal.dot(layerpars->
momentum().unit()));
2045 const double oldde = meff->
deltaE();
2047 std::unique_ptr<EnergyLoss> eloss;
2048 double sigmascat = 0;
2050 if (matprop !=
nullptr) {
2051 eloss = std::make_unique<EnergyLoss>(
2052 m_elosstool->energyLoss(*matprop, p, 1. / costracksurf,
2054 sigmascat = std::sqrt(
m_scattool->sigmaSquare(*matprop, p, 1. / costracksurf, matEffects));
2056 if (eloss !=
nullptr) {
2061 sigmascat = std::sqrt(
m_scattool->sigmaSquare(tmpprop, p, 1. / costracksurf, matEffects));
2065 sigmascat / std::sin(layerpars->parameters()[
Trk::theta]),
2079 }
else if (eloss !=
nullptr) {
2090 refpars[0], refpars[1], refpars[2], refpars[3], refpars[4], std::nullopt
2098 trajectory.
brems().clear();
2100 trajectory.
brems().resize(1);
2101 trajectory.
brems()[0] = bremdp;
2112 std::unique_ptr<Track> track(
myfit(ctx, cache, trajectory, *minpar, runOutlier, matEffects));
2114 bool pseudoupdated =
false;
2116 if ((track !=
nullptr) && hasid && hasmuon) {
2117 for (std::unique_ptr<GXFTrackState> & pseudostate : trajectory.
trackStates()) {
2119 (pseudostate ==
nullptr) ||
2121 pseudostate->fitQuality().chiSquared() < 10
2127 const std::unique_ptr<const TrackParameters> updpar(
m_updator->removeFromState(
2129 pseudostate->measurement()->localParameters(),
2130 pseudostate->measurement()->localCovariance()
2133 if (updpar ==
nullptr) {
2138 covMatrix(0, 0) = 100;
2140 std::unique_ptr<const PseudoMeasurementOnTrack> newpseudo = std::make_unique<const PseudoMeasurementOnTrack>(
2142 std::move(covMatrix),
2146 pseudostate->setMeasurement(std::move(newpseudo));
2148 errors[0] = errors[2] = errors[3] = errors[4] = -1;
2150 pseudostate->setMeasurementErrors(errors);
2151 pseudoupdated =
true;
2154 if (pseudoupdated) {
2157 track.reset(
myfit(ctx, cache, trajectory, *track->perigeeParameters(),
false,
muon));
2167 if (track !=
nullptr) {
2170 track->info().addPatternReco(old_info);
2173 return track.release();
2176 std::unique_ptr<Track>
2178 const EventContext& ctx,
2184 ATH_MSG_DEBUG(
"--> entering GlobalChi2Fitter::fit(PRDS,TP,)");
2187 for (
const auto *prd : prds) {
2188 const Surface & prdsurf = (*prd).detectorElement()->surface((*prd).identify());
2193 plsurf =
static_cast < const PlaneSurface *
>(&prdsurf);
2200 if ((slsurf ==
nullptr) && (plsurf ==
nullptr)) {
2201 ATH_MSG_ERROR(
"Surface is neither PlaneSurface nor StraightLineSurface!");
2206 }
else if (slsurf !=
nullptr) {
2215 }
else if (plsurf !=
nullptr) {
2216 if (param.covariance() !=
nullptr) {
2238 if (rot !=
nullptr) {
2239 rots.push_back(rot);
2243 std::unique_ptr<Track> track =
2244 fit(ctx, rots, param, runOutlier, matEffects);
2246 for (
const auto *rot : rots) {
2253 std::unique_ptr<Track>
2255 const EventContext& ctx,
2256 const Track& inputTrack,
2261 ATH_MSG_DEBUG(
"--> entering GlobalChi2Fitter::fit(Track,Meas'BaseSet,,)");
2276 if (minpar ==
nullptr) {
2294 for (; itStates != endState; ++itStates) {
2299 (trajectory.
trackStates().back()->materialEffects() !=
nullptr) &&
2300 trajectory.
trackStates().back()->materialEffects()->sigmaDeltaE() > 50.001
2302 trajectory.
trackStates().back()->materialEffects()->setKink(
true);
2308 for (
const auto & measBase : addMeasColl) {
2309 if (measBase ==
nullptr) {
2310 ATH_MSG_WARNING(
"There is an empty MeasurementBase object in the track! Skip this object..");
2318 std::unique_ptr<Track> track(
myfit(ctx, cache, trajectory, *minpar, runOutlier, matEffects));
2321 if (track !=
nullptr) {
2322 const double oldqual =
2327 const double newqual =
2328 track->fitQuality()->numberDoF() != 0 ?
2329 track->fitQuality()->chiSquared() / track->fitQuality()->numberDoF() :
2332 if (
m_extensioncuts && runOutlier && newqual > 2 && newqual > 2 * oldqual) {
2336 track.reset(
nullptr);
2340 if (track !=
nullptr) {
2349 std::unique_ptr<Track>
2351 const EventContext& ctx,
2357 ATH_MSG_DEBUG(
"--> entering GlobalChi2Fitter::fit(Track,PRDS,)");
2361 for (
const auto *prd : prds) {
2362 const Surface & prdsurf = (*prd).detectorElement()->surface((*prd).identify());
2365 std::unique_ptr<const TrackParameters>
const trackparForCorrect(
2381 rot =
m_ROTcreator->correct(*prd, *trackparForCorrect, ctx);
2384 if (rot !=
nullptr) {
2385 rots.push_back(rot);
2389 std::unique_ptr<Track> track =
fit(ctx,intrk, rots, runOutlier, matEffects);
2391 for (
const auto *rot : rots) {
2399 const EventContext& ctx,
2405 ATH_MSG_DEBUG(
"--> entering GlobalChi2Fitter::fit(Meas'BaseSet,,)");
2418 for (
const auto *itSet : rots) {
2419 if (itSet ==
nullptr) {
2420 ATH_MSG_WARNING(
"There is an empty MeasurementBase object in the track! Skip this object..");
2426 std::unique_ptr<const TrackParameters> startpar(param.
clone());
2429 matEffects ==
muon &&
2435 myfit(ctx,cache, trajectory, *startpar, runOutlier, matEffects);
2451 covMatrix(0, 0) = 100;
2453 std::unique_ptr<const PseudoMeasurementOnTrack> newpseudo = std::make_unique<const PseudoMeasurementOnTrack>(
2455 std::move(covMatrix),
2459 trajectory.
trackStates().front()->setMeasurement(std::move(newpseudo));
2464 covMatrix(0, 0) = 100;
2466 std::unique_ptr<const PseudoMeasurementOnTrack> newpseudo = std::make_unique<const PseudoMeasurementOnTrack>(
2468 std::move(covMatrix),
2472 trajectory.
trackStates().back()->setMeasurement(std::move(newpseudo));
2479 refpars[0], refpars[1], refpars[2], refpars[3], refpars[4], std::nullopt
2484 myfit(ctx,cache, trajectory, *startpar, runOutlier, matEffects);
2495 refpars[0], refpars[1], refpars[2], refpars[3], refpars[4], std::nullopt
2502 firstpar = trajectory.
trackStates().front()->trackParameters();
2505 covMatrix(0, 0) = 100;
2507 std::unique_ptr<const PseudoMeasurementOnTrack> newpseudo = std::make_unique<const PseudoMeasurementOnTrack>(
2509 std::move(covMatrix),
2513 trajectory.
trackStates().front()->setMeasurement(std::move(newpseudo));
2515 errors[0] = errors[2] = errors[3] = errors[4] = -1;
2517 trajectory.
trackStates().front()->setMeasurementErrors(errors);
2521 lastpar = trajectory.
trackStates().back()->trackParameters();
2524 covMatrix(0, 0) = 100;
2526 std::unique_ptr<const PseudoMeasurementOnTrack> newpseudo = std::make_unique<const PseudoMeasurementOnTrack>(
2528 std::move(covMatrix),
2532 trajectory.
trackStates().back()->setMeasurement(std::move(newpseudo));
2534 errors[0] = errors[2] = errors[3] = errors[4] = -1;
2536 trajectory.
trackStates().back()->setMeasurementErrors(errors);
2540 std::unique_ptr<Track> track;
2542 if (startpar !=
nullptr) {
2543 track.reset(
myfit(ctx,cache, trajectory, *startpar, runOutlier, matEffects));
2546 if (track !=
nullptr) {
2576 std::unique_ptr<GXFMaterialEffects> newmeff;
2584 newmeff = std::make_unique<GXFMaterialEffects>(*meff);
2588 const double sigmascat = std::sqrt(
m_scattool->sigmaSquare(
2609 newmeff = std::make_unique<GXFMaterialEffects>(newmeot);
2625 ((**trajectory.
trackStates().rbegin()).trackParameters() !=
nullptr)
2627 const double delta_p = 1000 * (
2629 (**trajectory.
trackStates().rbegin()).trackParameters()->
2633 newmeff->setdelta_p(delta_p);
2669 seg =
static_cast<const Segment *
>(measbase);
2679 for (
int i = 0; i <
imax; i++) {
2682 std::unique_ptr<GXFTrackState> ptsos = std::make_unique<GXFTrackState>(
2683 std::unique_ptr<const MeasurementBase>(measbase2->
clone()),
2684 std::unique_ptr<const TrackParameters>(newtrackpar !=
nullptr ? newtrackpar->
clone() :
nullptr)
2687 double sinstereo = 0;
2689 errors[0] = errors[1] = errors[2] = errors[3] = errors[4] = -1;
2700 bool measphi =
false;
2703 bool rotated =
false;
2706 if (
m_DetID->is_pixel(hitid)) {
2708 }
else if (
m_DetID->is_sct(hitid)) {
2709 if (covmat.cols() != 1 && covmat(1, 0) != 0) {
2713 }
else if (
m_DetID->is_trt(hitid)) {
2723 }
else if (
m_DetID->is_mdt(hitid)) {
2725 }
else if (
m_DetID->is_tgc(hitid)) {
2730 }
else if (
m_DetID->is_csc(hitid)) {
2732 }
else if (
m_DetID->is_mm(hitid)) {
2734 }
else if (
m_DetID->is_stgc(hitid)) {
2740 const auto [covEigenValueSmall, covStereoAngle] = principalComponentAnalysis2x2(covmat);
2741 errors[0] = std::sqrt(covEigenValueSmall);
2742 sinstereo = std::sin(covStereoAngle);
2744 errors[0] = std::sqrt(covmat(0, 0));
2746 errors[1] = std::sqrt(covmat(1, 1));
2757 const double dotprod1 = measdir.dot(
Amg::Vector3D(0, 0, 1));
2759 if (std::abs(dotprod1) < .5 && std::abs(dotprod2) < .5) {
2767 int param_index = 0;
2769 errors[0] = std::sqrt(covmat(0, 0));
2774 errors[1] = std::sqrt(covmat(param_index, param_index));
2779 errors[2] = std::sqrt(covmat(param_index, param_index));
2784 errors[3] = std::sqrt(covmat(param_index, param_index));
2789 errors[4] = std::sqrt(covmat(param_index, param_index));
2810 ptsos->setMeasurementErrors(errors);
2811 ptsos->setSinStereo(sinstereo);
2812 ptsos->setMeasurementType(hittype);
2813 ptsos->setMeasuresPhi(measphi);
2834 if (tvol ==
nullptr) {
2838 const Trk::BinnedArray < Trk::Layer > *confinedLayers = tvol->
confinedLayers();
2841 if (confinedLayers !=
nullptr) {
2843 for (
const auto & layer : confinedLayers->
arrayObjects()) {
2845 if (layer !=
nullptr) {
2850 if ((layIndex.
value() == 0) || (layer->layerMaterialProperties() ==
nullptr)) {
2861 disclay =
static_cast<const DiscLayer *
>(layer);
2864 if (disclay !=
nullptr) {
2865 if (disclay->
center().z() < 0) {
2870 }
else if (cyllay !=
nullptr) {
2881 for (
size_t ib = 0 ; ib < bsurf.size(); ++ib) {
2882 const Layer *layer = bsurf[ib]->surfaceRepresentation().materialLayer();
2884 if (layer ==
nullptr)
continue;
2888 if ((layIndex.
value() == 0) || (layer->layerMaterialProperties() ==
nullptr)) {
2895 cylsurf =
static_cast<const CylinderSurface *
>(&layer->surfaceRepresentation());
2900 discsurf =
static_cast<const DiscSurface *
>(&layer->surfaceRepresentation());
2902 if (discsurf !=
nullptr) {
2904 discsurf->
center().z() < 0 &&
2909 discsurf->
center().z() > 0 &&
2915 (cylsurf !=
nullptr) &&
2921 if ((cylsurf ==
nullptr) && (discsurf ==
nullptr)) {
2928 if (confinedVolumes !=
nullptr) {
2929 for (
const auto & volume : confinedVolumes->
arrayObjects()) {
2930 if (volume !=
nullptr) {
2946 std::optional<std::pair<Amg::Vector3D, double>>
2959 const double * pos = parforextrap.
position().data();
2962 const double sinphi = std::sin(parforextrap.parameters()[
Trk::phi0]);
2963 const double cosphi = std::cos(parforextrap.parameters()[
Trk::phi0]);
2964 const double sintheta = std::sin(parforextrap.parameters()[
Trk::theta]);
2965 const double costheta = std::cos(parforextrap.parameters()[
Trk::theta]);
2969 const double r = (std::abs(currentqoverp) > 1e-10) ? -sintheta / (currentqoverp * 300. * field[2]) : 1e6;
2970 const double xc = parforextrap.
position().x() -
r * sinphi;
2971 const double yc = parforextrap.
position().y() +
r * cosphi;
2972 const double phi0 = std::atan2(parforextrap.
position().y() - yc, parforextrap.
position().x() - xc);
2973 const double z0 = parforextrap.
position().z();
2974 const double delta_s = (surf.
center().
z() -
z0) / costheta;
2975 const double delta_phi = delta_s * sintheta /
r;
2976 const double x = xc + std::abs(
r) * std::cos(
phi0 + delta_phi);
2977 const double y = yc + std::abs(
r) * std::sin(
phi0 + delta_phi);
2979 const double perp = intersect.perp();
2986 const double costracksurf = std::abs(costheta);
2988 return std::make_pair(intersect, costracksurf);
2991 std::optional<std::pair<Amg::Vector3D, double>>
3006 const double * pos = parforextrap.
position().data();
3009 const double sinphi = std::sin(parforextrap.parameters()[
Trk::phi0]);
3010 const double cosphi = std::cos(parforextrap.parameters()[
Trk::phi0]);
3011 const double sintheta = std::sin(parforextrap.parameters()[
Trk::theta]);
3012 const double costheta = std::cos(parforextrap.parameters()[
Trk::theta]);
3013 const double tantheta = std::tan(parforextrap.parameters()[
Trk::theta]);
3014 const double r = (std::abs(currentqoverp) > 1e-10) ? -sintheta / (currentqoverp * 300. * field[2]) : 1e6;
3015 const double xc = parforextrap.
position().x() -
r * sinphi;
3016 const double yc = parforextrap.
position().y() +
r * cosphi;
3017 const double phi0 = std::atan2(parforextrap.
position().y() - yc, parforextrap.
position().x() - xc);
3018 const double z0 = parforextrap.
position().z();
3019 const double d = xc * xc + yc * yc;
3020 const double rcyl = surf.
bounds().
r();
3021 double mysqrt = ((
r + rcyl) * (
r + rcyl) - d) * (d - (
r - rcyl) * (
r - rcyl));
3027 mysqrt = std::sqrt(mysqrt);
3028 double firstterm = xc / 2 + (xc * (rcyl * rcyl -
r *
r)) / (2 * d);
3029 double secondterm = (mysqrt * yc) / (2 * d);
3030 const double x1 = firstterm + secondterm;
3031 const double x2 = firstterm - secondterm;
3032 firstterm = yc / 2 + (yc * (rcyl * rcyl -
r *
r)) / (2 * d);
3033 secondterm = (mysqrt * xc) / (2 * d);
3034 const double y1 = firstterm - secondterm;
3035 const double y2 = firstterm + secondterm;
3038 const double dist1 = (
x - x1) * (
x - x1) + (
y - y1) * (
y - y1);
3039 const double dist2 = (
x - x2) * (
x - x2) + (
y - y2) * (
y - y2);
3041 if (dist1 < dist2) {
3049 const double phi1 = std::atan2(
y - yc,
x - xc);
3052 const double delta_z =
r * deltaphi / tantheta;
3053 const double z =
z0 + delta_z;
3064 const Amg::Vector3D trackdir(cos(phidir) * sintheta, std::sin(phidir) * sintheta, costheta);
3066 const double costracksurf = std::abs(normal.unit().dot(trackdir));
3068 return std::make_pair(intersect, costracksurf);
3075 std::vector<std::pair<const Layer *, const Layer *>> &layers,
3080 std::vector<std::unique_ptr<GXFTrackState>> & states = trajectory.
trackStates();
3081 std::vector<std::unique_ptr<GXFTrackState>> oldstates = std::move(states);
3084 states.reserve(oldstates.size() + layers.size());
3092 for (
int i = 0; i <= indexoffset; i++) {
3102 for (
int i = indexoffset + 1; i < (int) oldstates.size(); i++) {
3103 const double rmeas = oldstates[i]->position().perp();
3104 const double zmeas = oldstates[i]->position().z();
3112 while (layerindex < (
int) layers.size()) {
3114 double costracksurf = 0.0;
3115 const Layer *layer =
nullptr;
3124 if (layers[layerindex].first !=
nullptr) {
3128 layer = layers[layerindex].first;
3135 if (oldstates[i]->trackParameters() !=
nullptr) {
3136 const double rlayer = cylsurf->
bounds().
r();
3137 if (std::abs(rmeas - rlayer) < std::abs(parforextrap->
position().perp() - rlayer)) {
3138 parforextrap = oldstates[i]->trackParameters();
3148 std::tie(intersect, costracksurf) =
res.value();
3154 if (cylsurf->
bounds().
r() > rmeas)
break;
3155 }
else if (layers[layerindex].second !=
nullptr) {
3161 layer = layers[layerindex].second;
3164 if (oldstates[i]->trackParameters() !=
nullptr) {
3165 const double zlayer = discsurf->
center().z();
3166 if (std::abs(zmeas - zlayer) < std::abs(parforextrap->
position().z() - zlayer)) {
3167 parforextrap = oldstates[i]->trackParameters();
3172 std::tie(intersect, costracksurf) =
res.value();
3178 if (std::abs(discsurf->
center().z()) > std::abs(zmeas))
break;
3180 throw std::logic_error(
"Unhandled surface.");
3187 const MaterialProperties *matprop = layer->layerMaterialProperties()->fullMaterial(intersect);
3188 if (matprop ==
nullptr) {
3201 const double actualx0 = X0 / costracksurf;
3202 const double de = -std::abs(
3206 (
m_p != 0.0 ? std::abs(
m_p) : std::abs(1. / currentqoverp)),
3208 const double sintheta = std::sin(parforextrap->parameters()[
Trk::theta]);
3209 const double sigmascat = std::sqrt(
m_scattool->sigmaSquare(
3211 (
m_p != 0.0 ? std::abs(
m_p) : std::abs(1. / currentqoverp)),
3215 std::unique_ptr<GXFMaterialEffects> meff = std::make_unique<GXFMaterialEffects>();
3216 meff->setDeltaE(de);
3217 meff->setScatteringSigmas(sigmascat / sintheta, sigmascat);
3218 meff->setX0(actualx0);
3219 meff->setSurface(&layer->surfaceRepresentation());
3220 meff->setMaterialProperties(matprop);
3226 std::unique_ptr<EnergyLoss> eloss;
3229 eloss = std::make_unique<EnergyLoss>(
m_elosstool->energyLoss(
3231 (
m_p != 0.0 ? std::abs(
m_p) : std::abs(1. / currentqoverp)),
3236 if (eloss !=
nullptr) {
3237 meff->setSigmaDeltaE(eloss->sigmaDeltaE());
3242 meff->setDeltaE(-5);
3244 meff->setScatteringSigmas(0, 0);
3247 meff->setSigmaDeltaE(50);
3248 if (eloss !=
nullptr) {
3249 meff->setSigmaDeltaEPos(eloss->sigmaPlusDeltaE());
3250 meff->setSigmaDeltaENeg(eloss->sigmaMinusDeltaE());
3255 "X0: " << meff->x0() <<
" qoverp: " << currentqoverp <<
3256 " sigmascat " << meff->sigmaDeltaTheta() <<
" eloss: " << meff->deltaE() <<
3257 " sigma eloss: " << meff->sigmaDeltaE()
3264 std::unique_ptr<GXFTrackState> matstate = std::make_unique<GXFTrackState>(
3266 std::unique_ptr<const TrackParameters>()
3268 matstate->setPosition(intersect);
3284 std::vector<std::pair<const Layer *, const Layer *>> & layers,
3285 std::vector<std::pair<const Layer *, const Layer *>> & upstreamlayers,
3286 const std::vector<std::unique_ptr<GXFTrackState>> & oldstates,
3295 upstreamlayers.reserve(5);
3314 const double lastz = laststate->
position().z();
3315 const double lastr = laststate->
position().perp();
3321 const double tantheta = std::tan(refpar->parameters()[
Trk::theta]);
3322 const double slope = (tantheta != 0) ? 1 / tantheta : 0;
3328 std::vector < const Layer *>::const_iterator it;
3329 std::vector < const Layer *>::const_iterator itend;
3347 for (; it != itend; ++it) {
3352 if (std::abs((*it)->surfaceRepresentation().center().z()) > std::abs(lastz)) {
3360 const DiscBounds *discbounds =
static_cast<const DiscBounds *
> (&(*it)->surfaceRepresentation().bounds());
3365 if (discbounds->
rMax() < firstr || discbounds->
rMin() > lastr) {
3369 const double rintersect = firstr + ((*it)->surfaceRepresentation().center().
z() - firstz) / slope;
3372 rintersect < discbounds->rMin() - 50 ||
3373 rintersect > discbounds->
rMax() + 50
3383 if ((*it) == endlayer) {
3395 std::abs((*it)->surfaceRepresentation().center().z()) < std::abs(firstz) ||
3398 upstreamlayers.emplace_back((
Layer *)
nullptr, (*it));
3405 (*it) != startlayer &&
3406 (std::abs((*it)->surfaceRepresentation().center().z()) > std::abs(firstz2) ||
3407 (*it) == startlayer2)
3409 layers.emplace_back((
Layer *)
nullptr, (*it));
3421 if (barrelcylinder->surfaceRepresentation().bounds().r() > lastr) {
3428 const double zintersect = firstz + (barrelcylinder->surfaceRepresentation().bounds().
r() - firstr) * slope;
3430 if (std::abs(zintersect - barrelcylinder->surfaceRepresentation().center().z()) >
3431 ((
const CylinderSurface*)(&barrelcylinder->surfaceRepresentation()))->bounds().halflengthZ() + 50) {
3435 if (barrelcylinder == endlayer) {
3442 if (barrelcylinder->surfaceRepresentation().bounds().r() < firstr ||
3443 barrelcylinder == startlayer) {
3444 upstreamlayers.emplace_back(barrelcylinder, (
Layer*)
nullptr);
3447 if (barrelcylinder != startlayer &&
3448 (barrelcylinder->surfaceRepresentation().bounds().r() > firstr2 ||
3449 barrelcylinder == startlayer2)) {
3450 layers.emplace_back(barrelcylinder, (
Layer*)
nullptr);
3463 const EventContext& ctx,
3474 if (!caloEntranceIsValid) {
3501 addMaterial(ctx, cache, trajectory, refpar2, matEffects);
3517 bool hasmat =
false;
3518 int indexoffset = 0;
3519 int lastmatindex = 0;
3520 std::vector<std::unique_ptr<GXFTrackState>> & oldstates = trajectory.
trackStates();
3534 for (
int i = 0; i < (int) oldstates.size(); i++) {
3535 if (oldstates[i]->materialEffects() !=
nullptr) {
3544 if (firstsistate ==
nullptr) {
3545 if (oldstates[i]->trackParameters() ==
nullptr) {
3546 std::unique_ptr<const TrackParameters> tmppar(
m_propagator->propagateParameters(
3549 oldstates[i]->associatedSurface(),
3556 if (tmppar ==
nullptr)
return;
3558 oldstates[i]->setTrackParameters(std::move(tmppar));
3560 firstsistate = oldstates[i].get();
3562 lastsistate = oldstates[i].get();
3570 if (lastsistate ==
nullptr) {
3571 throw std::logic_error(
"No track state");
3581 std::unique_ptr<const TrackParameters> tmppar(
m_propagator->propagateParameters(
3590 if (tmppar ==
nullptr)
return;
3602 indexoffset = lastmatindex;
3617 std::vector<std::pair<const Layer *, const Layer *>> layers;
3618 std::vector<std::pair<const Layer *, const Layer *>> & upstreamlayers = trajectory.
upstreamMaterialLayers();
3623 addMaterialGetLayers(cache, layers, upstreamlayers, oldstates, *firstsistate, *lastsistate, refpar, hasmat);
3632 const EventContext& ctx,
3638 if (refpar2 ==
nullptr) {
3648 std::vector<std::unique_ptr<GXFTrackState>> & states = trajectory.
trackStates();
3649 std::vector<std::unique_ptr<GXFTrackState>> matstates;
3650 std::unique_ptr< const std::vector < const TrackStateOnSurface *>,
3651 void (*)(
const std::vector<const TrackStateOnSurface *> *) >
3653 bool matvec_used=
false;
3654 std::unique_ptr<TrackParameters> startmatpar1;
3655 std::unique_ptr<TrackParameters> startmatpar2;
3666 int npseudomuon1 = 0;
3667 int npseudomuon2 = 0;
3669 for (
auto & state : states) {
3675 if (firstidhit ==
nullptr) {
3684 if (firsthit ==
nullptr) {
3685 firsthit = state->measurement();
3687 if (tp ==
nullptr) {
3697 if (tp ==
nullptr) {
3701 state->setTrackParameters(std::unique_ptr<const TrackParameters>(tp));
3707 lasthit = state->measurement();
3713 if (firstidhit ==
nullptr) {
3714 firstidhit = state->measurement();
3717 if ((firstidpar ==
nullptr) && (tp !=
nullptr)) {
3721 lastidhit = state->measurement();
3722 if (tp !=
nullptr) {
3727 if (firstsiliconpar ==
nullptr) {
3728 firstsiliconpar = tp;
3730 lastsiliconpar = tp;
3742 if (firstmuonhit ==
nullptr) {
3743 firstmuonhit = state->measurement();
3744 if (tp !=
nullptr) {
3748 lastmuonhit = state->measurement();
3749 if (tp !=
nullptr) {
3755 if (meff->
deltaE() == 0) {
3756 if (firstcalopar ==
nullptr) {
3757 firstcalopar = state->trackParameters();
3759 lastcalopar = state->trackParameters();
3761 if (firstmatpar ==
nullptr) {
3762 firstmatpar = state->trackParameters();
3767 std::unique_ptr<TrackParameters> refpar;
3768 AmgVector(5) newpars = refpar2->parameters();
3775 newpars[0], newpars[1], newpars[2], newpars[3], newpars[4], std::nullopt
3778 if (firstmatpar !=
nullptr) {
3783 if ((startmatpar1 ==
nullptr) || ((firstidhit !=
nullptr) && (firstmuonhit !=
nullptr))) {
3787 const double mass = trajectory.
mass();
3788 if (mass > 200 * MeV) {
3789 const AmgVector(5) & newpars = startmatpar2->parameters();
3790 const double oldp = std::abs(1 / newpars[
Trk::qOverP]);
3793 startmatpar2 = startmatpar2->associatedSurface().createUniqueTrackParameters(
3794 newpars[0], newpars[1], newpars[2], newpars[3],
3795 sign / std::sqrt(oldp * oldp + 2 * 100 * MeV * std::sqrt(oldp * oldp + mass * mass) + 100 * MeV * 100 * MeV),
3800 AmgVector(5) newpars = startmatpar1->parameters();
3803 startmatpar1 = startmatpar1->associatedSurface().createUniqueTrackParameters(
3804 newpars[0], newpars[1], newpars[2], newpars[3], newpars[4], std::nullopt
3807 newpars = startmatpar2->parameters();
3810 startmatpar2 = startmatpar2->associatedSurface().createUniqueTrackParameters(
3811 newpars[0], newpars[1], newpars[2], newpars[3], newpars[4], std::nullopt
3819 refpar->momentum().unit()
3822 const double distance = getDistance(distsol);
3825 ATH_MSG_DEBUG(
"Obtaining upstream layers from Extrapolator");
3828 std::unique_ptr<const TrackParameters> tmppar;
3830 if (firstmuonhit !=
nullptr) {
3832 if (caloEntranceIsValid) {
3839 if (tmppar !=
nullptr) {
3840 destsurf = &tmppar->associatedSurface();
3845 if (matvec_used) cache.
m_matTempStore.push_back( std::move(matvec) );
3850 false, matEffects) );
3853 if (matvec && !matvec->empty()) {
3854 for (
int i = (
int)matvec->size() - 1; i > -1; i--) {
3859 std::unique_ptr<GXFMaterialEffects> meff = std::make_unique<GXFMaterialEffects>(*meot);
3860 const TrackParameters * newpars = (*matvec)[i]->trackParameters() !=
nullptr ? (*matvec)[i]->trackParameters()->clone() :
nullptr;
3861 meff->setSigmaDeltaE(0);
3862 matstates.push_back(std::make_unique<GXFTrackState>(
3864 std::unique_ptr<const TrackParameters>(newpars)
3877 refpar->momentum().unit()
3880 const double distance = getDistance(distsol);
3883 ATH_MSG_DEBUG(
"Obtaining downstream ID layers from Extrapolator");
3885 std::unique_ptr<const TrackParameters> tmppar;
3886 std::unique_ptr<Surface> calosurf;
3887 if (firstmuonhit !=
nullptr) {
3889 if (caloEntranceIsValid) {
3897 if (tmppar !=
nullptr) {
3901 cylcalosurf =
static_cast<const CylinderSurface *
>(&tmppar->associatedSurface());
3906 disccalosurf =
static_cast<const DiscSurface *
>(&tmppar->associatedSurface());
3908 if (cylcalosurf !=
nullptr) {
3911 const double radius = cylbounds.
r();
3913 calosurf = std::make_unique<CylinderSurface>(trans, radius - 1, hlength);
3914 }
else if (disccalosurf !=
nullptr) {
3915 const double newz = (
3916 disccalosurf->
center().
z() > 0 ?
3917 disccalosurf->
center().z() - 1 :
3918 disccalosurf->
center().z() + 1
3922 disccalosurf->
center().x(),
3923 disccalosurf->
center().y(),
3928 trans.translation() << newpos;
3931 const double rmin = discbounds->
rMin();
3932 const double rmax = discbounds->
rMax();
3933 calosurf = std::make_unique<DiscSurface>(trans, rmin, rmax);
3935 destsurf = calosurf.release();
3939 if (matvec_used) cache.
m_matTempStore.push_back( std::move(matvec) );
3941 ctx, *startmatpar2, *destsurf,
alongMomentum,
false, matEffects));
3942 matvec_used =
false;
3944 if (matvec && !matvec->empty()) {
3945 for (
const auto & i : *matvec) {
3951 std::unique_ptr<GXFMaterialEffects> meff = std::make_unique<GXFMaterialEffects>(*meot);
3957 meff->setDeltaE(-5);
3960 meff->setScatteringSigmas(0, 0);
3963 meff->setSigmaDeltaE(50);
3966 const TrackParameters * newparams = i->trackParameters() !=
nullptr ? i->trackParameters()->clone() :
nullptr;
3968 matstates.push_back(std::make_unique<GXFTrackState>(
3970 std::unique_ptr<const TrackParameters>(newparams)
3982 if (cache.
m_calomat && (firstmuonhit !=
nullptr) && (firstidhit !=
nullptr)) {
3985 std::vector<MaterialEffectsOnTrack> calomeots =
m_calotool->extrapolationSurfacesAndEffects(
3994 if (calomeots.empty()) {
3999 std::unique_ptr<const TrackParameters> prevtrackpars =
unique_clone(lastidpar);
4000 if (lasthit == lastmuonhit) {
4001 for (
int i = 0; i < (int) calomeots.size(); i++) {
4004 std::unique_ptr<const TrackParameters> layerpar(
m_propagator->propagateParameters(
4007 calomeots[i].associatedSurface(),
4014 if (layerpar ==
nullptr) {
4019 std::unique_ptr<GXFMaterialEffects> meff = std::make_unique<GXFMaterialEffects>(calomeots[i]);
4022 lastcalopar = layerpar.get();
4026 const double qoverp = layerpar->parameters()[
Trk::qOverP];
4027 double qoverpbrem = 0;
4031 (firstmuonpar !=
nullptr) &&
4032 std::abs(firstmuonpar->parameters()[
Trk::qOverP]) > 1.e-9
4034 qoverpbrem = firstmuonpar->parameters()[
Trk::qOverP];
4036 const double sign = (qoverp > 0) ? 1 : -1;
4037 qoverpbrem =
sign / (1 / std::abs(qoverp) - std::abs(calomeots[i].energyLoss()->deltaE()));
4040 const AmgVector(5) & newpar = layerpar->parameters();
4042 layerpar = layerpar->associatedSurface().createUniqueTrackParameters(
4043 newpar[0], newpar[1], newpar[2], newpar[3], qoverpbrem, std::nullopt
4045 meff->setdelta_p(1000 * (qoverpbrem - qoverp));
4048 matstates.push_back(std::make_unique<GXFTrackState>(
4050 std::unique_ptr<const TrackParameters>(layerpar !=
nullptr ? layerpar->clone() :
nullptr)
4052 prevtrackpars = std::move(layerpar);
4057 firsthit == firstmuonhit &&
4061 for (
int i = 0; i < (int) calomeots.size(); i++) {
4063 std::unique_ptr<const TrackParameters> layerpar(
m_propagator->propagateParameters(
4066 calomeots[i].associatedSurface(),
4073 if (layerpar ==
nullptr) {
4078 std::unique_ptr<GXFMaterialEffects> meff = std::make_unique<GXFMaterialEffects>(calomeots[i]);
4087 const double qoverpbrem = layerpar->parameters()[
Trk::qOverP];
4092 (lastmuonpar !=
nullptr) &&
4093 std::abs(lastmuonpar->parameters()[
Trk::qOverP]) > 1.e-9
4097 const double sign = (qoverpbrem > 0) ? 1 : -1;
4098 qoverp =
sign / (1 / std::abs(qoverpbrem) + std::abs(calomeots[i].energyLoss()->deltaE()));
4101 meff->setdelta_p(1000 * (qoverpbrem - qoverp));
4102 const AmgVector(5) & newpar = layerpar->parameters();
4104 prevtrackpars = layerpar->associatedSurface().createUniqueTrackParameters(
4105 newpar[0], newpar[1], newpar[2], newpar[3], qoverp, std::nullopt
4109 matstates.insert(matstates.begin(), std::make_unique<GXFTrackState>(std::move(meff), std::move(layerpar)));
4114 if (lasthit == lastmuonhit && cache.
m_extmat) {
4115 std::unique_ptr<const Trk::TrackParameters> muonpar1;
4117 if (lastcalopar !=
nullptr) {
4119 if (msEntranceIsValid) {
4127 if (muonpar1 !=
nullptr) {
4135 rot.col(2) = trackdir;
4137 trans.linear().matrix() << rot;
4138 trans.translation() << muonpar1->position() - .1 * trackdir;
4141 std::unique_ptr<const TrackParameters> curvlinpar(
m_extrapolator->extrapolateDirectly(
4149 if (curvlinpar !=
nullptr) {
4150 muonpar1 = std::move(curvlinpar);
4154 muonpar1 = std::unique_ptr<TrackParameters>(lastcalopar->
clone());
4158 muonpar1 = std::unique_ptr<TrackParameters>(refpar->clone());
4163 if (muonpar1 !=
nullptr) {
4165 muonpar1->position(),
4166 muonpar1->momentum().unit()
4170 double distance = getDistance(distsol);
4173 0) and (firstmuonhit !=
nullptr)) {
4175 muonpar1->position(),
4176 muonpar1->momentum().unit()
4182 distance = distsol.
first();
4185 std::abs(distsol.
first()) < std::abs(distsol.
second()) ?
4191 if (distance < 0 && distsol.
numberOfSolutions() > 0 && (firstidhit ==
nullptr)) {
4192 if (firstmuonpar !=
nullptr) {
4193 AmgVector(5) newpars = firstmuonpar->parameters();
4200 newpars[0], newpars[1], newpars[2], newpars[3], newpars[4], std::nullopt
4203 std::unique_ptr<const TrackParameters> tmppar(
m_propagator->propagateParameters(
4213 if (tmppar !=
nullptr) {
4214 muonpar1 = std::move(tmppar);
4220 ATH_MSG_DEBUG(
"Obtaining downstream layers from Extrapolator");
4221 if (matvec_used) cache.
m_matTempStore.push_back( std::move(matvec) );
4224 states.back()->associatedSurface(),
4228 matvec_used =
false;
4235 if (matvec && !matvec->empty()) {
4236 for (
int j = 0; j < (int) matvec->size(); j++) {
4242 std::unique_ptr<GXFMaterialEffects> meff = std::make_unique<GXFMaterialEffects>(*meot);
4247 std::abs(meff->deltaE()) > 25 &&
4248 std::abs((*matvec)[j]->trackParameters()->position().z()) < 13000
4253 const TrackParameters * newparams = (*matvec)[j]->trackParameters() !=
nullptr ? (*matvec)[j]->trackParameters()->clone() :
nullptr;
4255 matstates.push_back(std::make_unique<GXFTrackState>(
4257 std::unique_ptr<const TrackParameters>(newparams)
4267 if (firsthit == firstmuonhit && cache.
m_extmat && (firstcalopar !=
nullptr)) {
4268 std::unique_ptr<const Trk::TrackParameters> muonpar1;
4271 if (msEntranceIsValid) {
4279 if (muonpar1 !=
nullptr) {
4287 rot.col(2) = trackdir;
4289 trans.linear().matrix() << rot;
4290 trans.translation() << muonpar1->position() - .1 * trackdir;
4293 std::unique_ptr<const TrackParameters> curvlinpar(
m_extrapolator->extrapolateDirectly(
4301 if (curvlinpar !=
nullptr) {
4302 muonpar1 = std::move(curvlinpar);
4306 muonpar1 = std::unique_ptr<const TrackParameters>(firstcalopar->
clone());
4312 if (muonpar1 !=
nullptr) {
4314 muonpar1->position(),
4315 muonpar1->momentum().unit()
4319 const double distance = getDistance(distsol);
4323 ATH_MSG_DEBUG(
"Collecting upstream muon material from extrapolator");
4324 if (matvec_used) cache.
m_matTempStore.push_back( std::move(matvec) );
4327 states[0]->associatedSurface(),
4331 matvec_used =
false;
4333 if (matvec && !matvec->empty()) {
4334 ATH_MSG_DEBUG(
"Retrieved " << matvec->size() <<
" material states");
4336 for (
int j = 0; j < (int) matvec->size(); j++) {
4339 if (meb !=
nullptr) {
4344 std::unique_ptr<GXFMaterialEffects> meff = std::make_unique<GXFMaterialEffects>(*meot);
4349 std::abs(meff->deltaE()) > 25 &&
4350 std::abs((*matvec)[j]->trackParameters()->position().z()) < 13000
4356 (*matvec)[j]->trackParameters() !=
nullptr
4357 ? (*matvec)[j]->trackParameters()->clone()
4360 matstates.insert(matstates.begin(), std::make_unique<GXFTrackState>(
4362 std::unique_ptr<const TrackParameters>(tmpparams)
4375 std::vector<std::unique_ptr<GXFTrackState>> & newstates = states;
4376 std::vector<std::unique_ptr<GXFTrackState>> oldstates = std::move(newstates);
4379 newstates.reserve(oldstates.size() + matstates.size());
4382 int firstlayerno = -1;
4388 const double cosphi = std::cos(refpar->parameters()[
Trk::phi0]);
4389 const double sinphi = std::sin(refpar->parameters()[
Trk::phi0]);
4391 for (
int i = cache.
m_acceleration ? 1 : 0; i < (
int) oldstates.size(); i++) {
4392 bool addlayer =
true;
4394 while (addlayer && layerno < (
int) matstates.size()) {
4396 const TrackParameters *layerpar = matstates[layerno]->trackParameters();
4398 const DistanceSolution distsol = oldstates[i]->associatedSurface().straightLineDistanceEstimate(
4403 const double distance = getDistance(distsol);
4411 const double trackimpact = std::abs(-refpar->position().x() * sinphi + refpar->position().y() * cosphi);
4413 if (trackimpact > cylinderradius - 5 * mm) {
4419 if (i == (
int) oldstates.size() - 1) {
4428 ATH_MSG_DEBUG(
"Material state with excessive scattering, skipping it");
4439 if (firstlayerno < 0) {
4440 firstlayerno = layerno;
4445 if ((layerpar !=
nullptr) && matEffects !=
pion && matEffects !=
muon) {
4447 const Layer *lay =
nullptr;
4449 if (tvol !=
nullptr) {
4465 if (matvec_used) cache.
m_matTempStore.push_back( std::move(matvec) );
4479 const AmgVector(5) & pars = param.parameters();
4481 pars[0], pars[1], pars[2], pars[3], pars[4], std::nullopt
4490 std::unique_ptr<const TrackParameters> per(
m_extrapolator->extrapolate(
4491 Gaudi::Hive::currentContext(),param, tmppersf,
oppositeMomentum,
false, matEffects));
4493 if (per ==
nullptr) {
4494 ATH_MSG_DEBUG(
"Cannot make Perigee with starting parameters");
4498 if(std::abs(per->position().z())>5000.) {
4499 ATH_MSG_WARNING(
"Pathological perigee well outside of tracking detector!! Returning nullptr");
4507 const EventContext& ctx,
4514 ATH_MSG_DEBUG(
"--> entering GlobalChi2Fitter::myfit_helper");
4542 if (trajectory.
nDOF() < 0) {
4552 ATH_MSG_WARNING(
"Attempt to apply material corrections with q/p=0, reject track");
4566 std::unique_ptr<const TrackParameters> per =
makePerigee(cache, param, matEffects);
4582 ATH_MSG_WARNING(
"Tracking Geometry Service and/or Material Updator Tool not configured");
4593 addMaterial(ctx, cache, trajectory, per !=
nullptr ? per.get() : ¶m, matEffects);
4596 ctx, cache, trajectory, per !=
nullptr ? per.get() : ¶m, matEffects);
4608 for (
const auto & state : trajectory.
trackStates()) {
4612 state->materialEffects()->deltaE() == 0
4614 scatstate2 = state.get();
4619 scatstate = state.get();
4626 if ((scatstate ==
nullptr) || (scatstate2 ==
nullptr)) {
4627 throw std::logic_error(
"Invalid scatterer");
4632 const int nstates = (int) trajectory.
trackStates().size();
4634 trajectory.
trackStates()[nstates / 2 - 1]->position() +
4640 std::unique_ptr<const TrackParameters> nearestpar;
4641 double mindist = 99999;
4642 std::vector < GXFTrackState * >mymatvec;
4645 if ((*it).trackParameters() ==
nullptr) {
4649 const double distance = persurf
4651 (*it).trackParameters()->position(),
4652 (*it).trackParameters()->momentum().unit())
4655 const bool insideid = (
4661 (((*it).measurement() !=
nullptr) && insideid) || (
4662 ((*it).materialEffects() !=
nullptr) &&
4664 (*it).materialEffects()->deltaE() == 0 ||
4665 ((*it).materialEffects()->sigmaDeltaPhi() == 0 &&
4667 (*it).materialEffects()->deltaPhi() != 0
4671 const double dist = ((*it).trackParameters()->position() -
vertex).perp();
4672 if (dist < mindist) {
4680 if (((*it).materialEffects() !=
nullptr) && distance > 0) {
4681 mymatvec.push_back(it.get());
4685 if (nearestpar ==
nullptr) {
4689 for (
auto & state : mymatvec) {
4691 const Surface &matsurf = state->associatedSurface();
4693 nearestpar->position(), nearestpar->momentum().unit());
4695 const double distance = getDistance(distsol);
4701 std::unique_ptr<const TrackParameters> tmppar(
m_propagator->propagateParameters(
4711 if (tmppar ==
nullptr) {
4723 if (tmppar ==
nullptr) {
4733 AmgVector(5) newpars = tmppar->parameters();
4735 if (state->materialEffects()->sigmaDeltaE() > 0) {
4736 newpars[
Trk::qOverP] += .001 * state->materialEffects()->delta_p();
4739 const double de = std::abs(state->materialEffects()->deltaE());
4740 const double oldp = std::abs(1 / newpars[
Trk::qOverP]);
4741 const double newp2 = oldp * oldp - 2 * de * std::sqrt(mass * mass + oldp * oldp) + de * de;
4747 nearestpar = tmppar->associatedSurface().createUniqueTrackParameters(
4748 newpars[0], newpars[1], newpars[2], newpars[3], newpars[4], std::nullopt
4752 std::unique_ptr<Trk::TrackParameters> tmpPars(
m_propagator->propagateParameters(
4763 if (tmpPars !=
nullptr) {
4764 per.reset(
static_cast < const Perigee *
>(tmpPars.release()));
4769 const double oldp = 1. / std::abs(per->parameters()[
Trk::qOverP]);
4771 const double newp = std::sqrt(oldp * oldp + 2 * toteloss * std::sqrt(oldp * oldp + mass * mass) + toteloss * toteloss);
4772 AmgVector(5) params = per->parameters();
4775 per = per->associatedSurface().createUniqueTrackParameters(
4776 params[0], params[1], params[2], params[3], params[4], std::nullopt
4780 if (per ==
nullptr) {
4797 }
else if (per ==
nullptr) {
4814 const AmgVector(5) & pars = per->parameters();
4815 per = per->associatedSurface().createUniqueTrackParameters(
4816 pars[0], pars[1], pars[2], pars[3], 0, std::nullopt
4822 if (per !=
nullptr) {
4832 Eigen::MatrixXd a_inv;
4833 a.resize(nfitpar, nfitpar);
4838 derivPool.setZero();
4840 for (std::unique_ptr<GXFTrackState> & state : trajectory.
trackStates()) {
4841 if (state->materialEffects() !=
nullptr) {
4844 state->setDerivatives(derivPool);
4847 bool doderiv =
true;
4850 for (
int it = 0; it <
m_maxit; ++it) {
4863 runIteration(ctx, cache, trajectory, it,
a, b, lu, doderiv);
4879 const double redchi2 = (trajectory.
nDOF() > 0) ? trajectory.
chi2() / trajectory.
nDOF() : 0;
4880 const double prevredchi2 = (trajectory.
nDOF() > 0) ? trajectory.
prevchi2() / trajectory.
nDOF() : 0;
4889 (redchi2 < prevredchi2 &&
4890 (redchi2 > prevredchi2 - 1 || redchi2 < 2)) ||
4891 nsihits + ntrthits == nhits
4896 if (it != 1 || nsihits != 0 || trajectory.
nDOF() <= 0 || trajectory.
chi2() / trajectory.
nDOF() <= 3) {
4911 if (ntrtprechits+ntrttubehits) {
4912 phf = float(ntrtprechits)/float(ntrtprechits+ntrttubehits);
4914 if (phf<m_minphfcut && it>=3) {
4915 if ((ntrtprechits+ntrttubehits)>=15) {
4919 ATH_MSG_DEBUG(
"Iter = " << it <<
" | nTRTStates = " << ntrthits
4920 <<
" | nTRTPrecHits = " << ntrtprechits
4921 <<
" | nTRTTubeHits = " << ntrttubehits
4922 <<
" | nOutliers = "
4942 if (trajectory.
prefit() == 0) {
4946 Eigen::LLT < Eigen::MatrixXd >
const lltOfW(
a);
4947 if (lltOfW.info() == Eigen::Success) {
4951 const int ncols =
a.cols();
4952 Amg::MatrixX const weightInvAMG = Amg::MatrixX::Identity(ncols, ncols);
4953 a_inv = lltOfW.solve(weightInvAMG);
4972 if (traj != &trajectory) {
4977 finaltrajectory = traj;
4985 if (nperpars == 5) {
4986 for (
int i = 0; i <
a.cols(); i++) {
4987 a_inv(4, i) *= .001;
4988 a_inv(i, 4) *= .001;
4992 int scatterPos = nperpars + 2 * nscat;
4993 for (
int bremno = 0; bremno < nbrem; bremno++, scatterPos++) {
4994 for (
int i = 0; i <
a.cols(); i++) {
4995 a_inv(scatterPos, i) *= .001;
4996 a_inv(i, scatterPos) *= .001;
5003 for (
int i = 0; i < nperparams; i++) {
5004 for (
int j = 0; j < nperparams; j++) {
5005 (errmat) (j, i) = a_inv(j, i);
5010 (errmat) (4, 4) = 1e-20;
5014 std::unique_ptr<const TrackParameters> measper(
5016 perpars[0], perpars[1], perpars[2], perpars[3], perpars[4], std::move(errmat)
5026 std::unique_ptr<Track> track =
nullptr;
5028 if (finaltrajectory->
prefit() > 0) {
5029 if (finaltrajectory != &trajectory) {
5031 delete finaltrajectory;
5037 track =
makeTrack(ctx,cache, *finaltrajectory, matEffects);
5050 (track !=
nullptr) && (
5051 track->fitQuality()->numberDoF() != 0 &&
5052 track->fitQuality()->chiSquared() / track->fitQuality()->numberDoF() > cut
5055 track.reset(
nullptr);
5059 if (track ==
nullptr) {
5063 if (finaltrajectory != &trajectory) {
5064 delete finaltrajectory;
5067 return track.release();
5071 const EventContext& ctx,
5072 const Cache & cache,
5076 int & bremno_maxbrempull,
5081 std::vector<std::unique_ptr<GXFTrackState>> & states = trajectory.
trackStates();
5101 const int nmeas = (int)
res.size();
5111 const int nDOF = trajectory.
nDOF();
5112 const bool doNewPseudoMeasurements = (
5116 std::abs((trajectory.
prevchi2() - trajectory.
chi2()) / nDOF) < 15 &&
5118 (nperpars == 0 || nidhits > 0)
5129 double maxbrempull = -0.2;
5138 for (
int hitno = 0; hitno < (int) states.size(); hitno++) {
5139 std::unique_ptr<GXFTrackState> & state = states[hitno];
5153 doNewPseudoMeasurements &&
5155 !state->associatedSurface().isFree() &&
5156 !state->isRecalibrated()
5159 covMatrix(0, 0) = 100;
5161 std::unique_ptr<const PseudoMeasurementOnTrack> newpseudo = std::make_unique<const PseudoMeasurementOnTrack>(
5163 std::move(covMatrix),
5167 state->setMeasurement(std::move(newpseudo));
5168 measbase = state->measurement();
5175 double *errors = state->measurementErrors();
5177 for (
int i = 0; i < 5; i++) {
5197 res[measno] = residuals[i];
5203 res[measno] = -std::remainder(-
res[measno], 2 *
M_PI);
5212 double *errors = state->measurementErrors();
5213 for (
int i = 0; i < 5; i++) {
5214 if (errors[i] > 0) {
5215 error[measno] = errors[i];
5226 ((trajectory.
prefit() == 0) || state->materialEffects()->deltaE() == 0)
5228 const double deltaPhi = state->materialEffects()->deltaPhi();
5229 const double measDeltaPhi = state->materialEffects()->measuredDeltaPhi();
5230 const double sigma2deltaPhi = std::pow(state->materialEffects()->sigmaDeltaPhi(), 2);
5231 const double deltaTheta = state->materialEffects()->deltaTheta();
5232 const double sigma2deltaTheta = std::pow(state->materialEffects()->sigmaDeltaTheta(), 2);
5234 if (trajectory.
prefit() != 1) {
5235 b[nperpars + 2 * scatno] -= (
deltaPhi - measDeltaPhi) / sigma2deltaPhi;
5236 b[nperpars + 2 * scatno + 1] -= deltaTheta / sigma2deltaTheta;
5238 b[nperpars + scatno] -= deltaTheta / sigma2deltaTheta;
5243 deltaTheta * deltaTheta / sigma2deltaTheta
5252 if ((state->materialEffects() !=
nullptr) && state->materialEffects()->sigmaDeltaE() > 0) {
5253 double averagenergyloss = std::abs(state->materialEffects()->deltaE());
5254 const double qoverpbrem = limitInversePValue(1000 * states[hitno]->trackParameters()->parameters()[
Trk::qOverP]);
5255 const double qoverp = limitInversePValue(qoverpbrem - state->materialEffects()->delta_p());
5256 const double pbrem = 1. / std::abs(qoverpbrem);
5257 const double p = 1. / std::abs(qoverp);
5258 const double mass = .001 * trajectory.
mass();
5259 const double energy = std::sqrt(p * p + mass * mass);
5260 const double bremEnergy = std::sqrt(pbrem * pbrem + mass * mass);
5262 const double resMaterial = .001 * averagenergyloss - energy + bremEnergy;
5263 res[nmeas - nbrem + bremno] = resMaterial;
5265 const double sigde = state->materialEffects()->sigmaDeltaE();
5266 const double sigdepos = state->materialEffects()->sigmaDeltaEPos();
5267 const double sigdeneg = state->materialEffects()->sigmaDeltaENeg();
5269 double errorMaterial = .001 * state->materialEffects()->sigmaDeltaE();
5270 error[nmeas - nbrem + bremno] = errorMaterial;
5281 if (state->materialEffects()->isKink()) {
5282 maxbrempull = -999999999;
5283 state_maxbrempull =
nullptr;
5289 trajectory.
prefit() == 0 &&
5291 sigde != sigdepos &&
5294 const double elosspull = resMaterial / errorMaterial;
5296 if (trajectory.
mass() > 100) {
5301 if (std::abs(elosspull) > 1) {
5302 if (elosspull < -1) {
5303 state->materialEffects()->setSigmaDeltaE(sigdepos);
5305 state->materialEffects()->setSigmaDeltaE(sigdeneg);
5308 errorMaterial = .001 * state->materialEffects()->sigmaDeltaE();
5309 error[nmeas - nbrem + bremno] = errorMaterial;
5321 !state->materialEffects()->isKink() && (
5322 (
m_fixbrem == -1 && elosspull < maxbrempull) ||
5326 bremno_maxbrempull = bremno;
5327 state_maxbrempull = state.get();
5328 maxbrempull = elosspull;
5337 trajectory.
prefit() == 0 &&
5338 state->materialEffects()->sigmaDeltaPhi() == 0 &&
5339 state->materialEffects()->isMeasuredEloss() &&
5340 resMaterial / (.001 * state->materialEffects()->sigmaDeltaEAve()) > 2.5
5342 const TrackParameters* parforcalo = states[hitno - 2]->trackParameters();
5345 std::vector<MaterialEffectsOnTrack> calomeots =
5358 if (calomeots.size() == 3) {
5359 averagenergyloss = std::abs(calomeots[1].energyLoss()->deltaE());
5360 const double newres = .001 * averagenergyloss - energy + bremEnergy;
5361 const double newerr = .001 * calomeots[1].energyLoss()->sigmaDeltaE();
5363 const double oldPull = resMaterial / errorMaterial;
5364 const double newPull = newres / newerr;
5366 if (std::abs(newPull) < std::abs(oldPull)) {
5367 ATH_MSG_DEBUG(
"Changing from measured to parametrized energy loss");
5369 state->materialEffects()->setEloss(std::unique_ptr<EnergyLoss>(calomeots[1].energyLoss()->clone()));
5370 state->materialEffects()->setSigmaDeltaE(calomeots[1].energyLoss()->sigmaDeltaE());
5371 res[nmeas - nbrem + bremno] = newres;
5372 error[nmeas - nbrem + bremno] = newerr;
5376 state->materialEffects()->setMeasuredEloss(
false);
5386 for (
int imeas = 0; imeas < nmeas; imeas++) {
5387 if (
error[imeas] == 0) {
5402 const Cache & cache,
5408 const double oldChi2 = trajectory.
prevchi2();
5409 const double newChi2 = trajectory.
chi2();
5414 const double nDOF = trajectory.
nDOF();
5415 const double oldRedChi2 = (nDOF > 0) ? oldChi2 / nDOF : 0;
5416 const double newRedChi2 = (nDOF > 0) ? newChi2 / nDOF : 0;
5419 trajectory.
prefit() > 0 && (
5420 (newRedChi2 < 2 && it != 0) ||
5421 (newRedChi2 < oldRedChi2 + .1 && std::abs(newRedChi2 - oldRedChi2) < 1 && it != 1)
5434 int miniter = (nsihits != 0 && nsihits + ntrthits == nhits) ? 1 : 2;
5435 miniter = std::max(miniter, cache.
m_miniter);
5437 if (it >= miniter && std::abs(oldChi2 - newChi2) < 1) {
5444 const int bremno_maxbrempull,
5450 if (state_maxbrempull ==
nullptr) {
5462 const int nmeas = (int)
res.size();
5465 const double oldError =
error[nmeas - nbrem + bremno_maxbrempull];
5467 error[nmeas - nbrem + bremno_maxbrempull] = newError;
5470 if (
a.cols() != nFitPars) {
5474 const double errorRatio = oldError / newError;
5475 const double errorReductionRatio = 1 - std::pow(errorRatio, 2);
5478 for (
int i = 0; i < nFitPars; i++) {
5479 if (weightderiv(nmeas - nbrem + bremno_maxbrempull, i) == 0) {
5483 for (
int j = i; j < nFitPars; j++) {
5484 const double newaij =
a(i, j) - errorReductionRatio *
5485 weightderiv(nmeas - nbrem + bremno_maxbrempull, i) *
5486 weightderiv(nmeas - nbrem + bremno_maxbrempull, j);
5488 a.fillSymmetric(i, j, newaij);
5490 weightderiv(nmeas - nbrem + bremno_maxbrempull, i) *= errorRatio;
5499 std::vector<std::unique_ptr<GXFTrackState>> & states = trajectory.
trackStates();
5512 const int nmeas = (int) weightderiv.rows();
5514 for (std::unique_ptr<GXFTrackState> & state : states) {
5518 const auto [scatmin, scatmax] = std::minmax(scatno, nscatupstream);
5519 const auto [bremmin, bremmax] = std::minmax(bremno, nbremupstream);
5526 const double sinStereo =
5528 state->sinStereo() :
5530 const double cosStereo =
5532 std::sqrt(1 - std::pow(sinStereo, 2)) :
5540 auto getThisDeriv = [sinStereo, cosStereo, &derivatives](
int i,
int j) ->
double {
5541 if (i == 0 && sinStereo != 0) {
5542 return derivatives(0, j) * cosStereo + sinStereo * derivatives(1, j);
5544 return derivatives(i, j);
5548 for (
int i = 0; i < 5; i++) {
5563 if (i == 0 && sinStereo != 0) {
5564 weightderiv.row(measno).head(cols) =
5565 (derivatives.row(0).
head(cols) * cosStereo +
5566 sinStereo * derivatives.row(1).head(cols)) /
5569 weightderiv.row(measno).head(cols) = derivatives.row(i).head(cols) /
error[measno];
5573 for (
int j = scatmin; j < scatmax; j++) {
5574 if (trajectory.
prefit() == 1) {
5575 const int index = nperparams + j;
5578 const int index = nperparams + 2 * j;
5580 weightderiv(measno,
index + 1) = getThisDeriv(i,
index + 1) /
error[measno];
5584 for (
int j = bremmin; j < bremmax; j++) {
5585 const int index = j + nperparams + 2 * nscat;
5592 double *errors = state->measurementErrors();
5593 for (
int i = 0; i < 5; i++) {
5594 if (errors[i] > 0) {
5600 ((trajectory.
prefit() == 0) || state->materialEffects()->deltaE() == 0)
5605 if ((state->materialEffects() !=
nullptr) && state->materialEffects()->sigmaDeltaE() > 0) {
5607 const double qoverpbrem = limitInversePValue(1000 * state->trackParameters()->parameters()[
Trk::qOverP]);
5608 const double qoverp = limitInversePValue(qoverpbrem - state->materialEffects()->delta_p());
5610 const double mass = .001 * trajectory.
mass();
5612 const auto thisMeasurementIdx{nmeas - nbrem + bremno};
5616 auto multiplier = [] (
double mass,
double qOverP){
5619 const auto qoverpTerm {multiplier(mass, qoverp) /
error[thisMeasurementIdx]};
5620 const auto qoverpBremTerm {multiplier(mass, qoverpbrem) /
error[thisMeasurementIdx]};
5623 weightderiv(thisMeasurementIdx, 4) = qoverpBremTerm - qoverpTerm;
5626 const auto bremNoBase = nperparams + 2 * nscat;
5627 if (bremno < nbremupstream) {
5628 weightderiv(thisMeasurementIdx, bremNoBase + bremno) = qoverpTerm;
5629 for (
int bremno2 = bremno + 1; bremno2 < nbremupstream; bremno2++) {
5630 weightderiv(thisMeasurementIdx, bremNoBase + bremno2) = qoverpTerm - qoverpBremTerm;
5633 weightderiv(thisMeasurementIdx, bremNoBase + bremno) = qoverpBremTerm;
5634 for (
int bremno2 = nbremupstream; bremno2 < bremno; bremno2++) {
5635 weightderiv(thisMeasurementIdx, bremNoBase + bremno2) = qoverpBremTerm - qoverpTerm;
5654 const int nMeas = (int)
res.size();
5659 for (
int i = 0; i < nPerPars; i++) {
5667 for (
int i = 0; i < (int) trajectory.
trackStates().size(); i++) {
5668 const std::unique_ptr<GXFTrackState> & state = trajectory.
trackStates()[i];
5671 if (meff ==
nullptr) {
5672 measno += state->numberOfMeasuredParameters();
5676 const int firstMeasurement = i < nUpstreamStates ? 0 : measno;
5677 const int lastMeasurement = i < nUpstreamStates ? measno : nMeas - nBrem;
5680 && (trajectory.
prefit() == 0 || meff->
deltaE() == 0)) {
5681 const int scatterPos = nPerPars + 2 * scatno;
5693 const int bremPos = nPerPars + nScatPars + bremno;
5704 const Cache & cache,
5717 const int nMeas = (int)
res.size();
5719 for (
int k = 0; k < nFitPars; k++) {
5728 for (
int measno = minMeasK; measno < maxMeasK; measno++) {
5729 b[k] +=
res[measno] * (1. /
error[measno]) * weightDeriv(measno, k);
5739 if (k == 4 || k >= nPerPars + nScatPars) {
5740 for (
int measno = nMeas - nBrem; measno < nMeas; measno++) {
5741 b[k] +=
res[measno] * (1. /
error[measno]) * weightDeriv(measno, k);
5748 const Cache & cache,
5755 for (
int k = 0; k < nFitPars; k++) {
5756 for (
int l = k; l < nFitPars; l++) {
5761 for (
int measno = minMeas; measno < maxMeas; measno++) {
5762 a_kl += weightDeriv(measno, k) * weightDeriv(measno, l);
5765 a.fillSymmetric(l, k, a_kl);
5783 const int nMeas = (int)
res.size();
5790 for (
int k = nPerPars; k < nPerPars + nScatPars; k += 2) {
5791 a(k, k) += 1. / std::pow(scatSigmas[scatno].first, 2);
5792 a(k + 1, k + 1) += 1. / std::pow(scatSigmas[scatno].second, 2);
5800 for (
int measno = nMeas - nBrem; measno < nMeas; measno++) {
5801 for (
int k = 4; k < nFitPars; k++) {
5803 k = nPerPars + nScatPars;
5806 for (
int l = k; l < nFitPars; l++) {
5808 l = nPerPars + nScatPars;
5811 const double a_kl =
a(l, k) + weightDeriv(measno, k) * weightDeriv(measno, l);
5812 a.fillSymmetric(l, k, a_kl);
5824 const double oldRedChi2,
5825 const double newRedChi2
5833 bool weightChanged =
false;
5840 double newPhiWeight = 1.1;
5841 double newThetaWeight = 1.001;
5842 if (trajectory.
prefit() == 0) {
5848 newPhiWeight = 1.00000001;
5849 }
else if (it == 1) {
5850 newPhiWeight = 1.0000001;
5851 }
else if (it <= 3) {
5852 newPhiWeight = 1.0001;
5853 }
else if (it <= 6) {
5854 newPhiWeight = 1.01;
5857 if (newRedChi2 > oldRedChi2 - 1 && newRedChi2 < oldRedChi2) {
5858 newPhiWeight = 1.0001;
5859 newThetaWeight = 1.0001;
5860 }
else if (newRedChi2 > oldRedChi2 - 25 && newRedChi2 < oldRedChi2) {
5861 newPhiWeight = 1.001;
5862 newThetaWeight = 1.0001;
5869 std::size_t scatno = 0;
5874 for (
const auto & state : trajectory.
trackStates()) {
5877 if (meff ==
nullptr) {
5881 const bool isValidPlaneSurface =
5883 static_cast<const PlaneSurface *
>(&state->associatedSurface()) !=
nullptr;
5888 if (meff->
deltaE() == 0 || (trajectory.
prefit() == 0 && isValidPlaneSurface)) {
5889 weightChanged =
true;
5891 const int scatNoIndex = 2 * scatno + nPerPars;
5895 std::stringstream message;
5896 message <<
"scatno is out of range " << scatno <<
" !< " << cache.
m_phiweight.size();
5897 throw std::range_error(message.str());
5905 a(scatNoIndex, scatNoIndex) /= cache.
m_phiweight[scatno];
5909 a(scatNoIndex, scatNoIndex) *= newPhiWeight;
5910 }
else if (trajectory.
prefit() >= 2) {
5911 a(scatNoIndex, scatNoIndex) *= newPhiWeight;
5912 a(scatNoIndex + 1, scatNoIndex + 1) *= newThetaWeight;
5939 trajectory.
prefit() == 2 &&
5942 (newRedChi2 < oldRedChi2 - 25 || newRedChi2 > oldRedChi2)
5947 return weightChanged;
5956 std::size_t scatno = 0;
5966 std::stringstream message;
5967 message <<
"scatno is out of range " << scatno <<
" !< " << cache.
m_phiweight.size();
5968 throw std::range_error(message.str());
5971 const bool isValidPlaneSurface =
5973 static_cast<const PlaneSurface *
>(&state->associatedSurface()) !=
nullptr;
5975 if (meff->
deltaE() == 0 || isValidPlaneSurface) {
5976 const int scatNoIndex = 2 * scatno + nPerPars;
5977 a(scatNoIndex, scatNoIndex) /= cache.
m_phiweight[scatno];
5992 const EventContext& ctx,
6001 const int nDOFold = trajectory.
nDOF();
6002 const double oldChi2 = trajectory.
chi2();
6003 const double oldRedChi2 = nDOFold > 0 ? oldChi2 / nDOFold : 0;
6025 int bremno_maxbrempull = 0;
6041 if ((state_maxbrempull !=
nullptr) && trajectory.
converged()) {
6050 const int nDOFnew = trajectory.
nDOF();
6051 const double newChi2 = trajectory.
chi2();
6052 const double newRedChi2 = nDOFnew > 0 ? newChi2 / nDOFnew : 0;
6054 ATH_MSG_DEBUG(
"old chi2: " << oldChi2 <<
"/" << nDOFold <<
"=" << oldRedChi2 <<
6055 ", new chi2: " << newChi2 <<
"/" << nDOFnew <<
"=" << newRedChi2);
6090 if (doDeriv || weightChanged) {
6101 if (trajectory.
prefit() == 0) {
6107 if (nSiHits + nTrtHits !=
nHits) {
6114 (newRedChi2 < 2 || (newRedChi2 < oldRedChi2 && newRedChi2 > oldRedChi2 - .5))
6136 Eigen::LLT<Eigen::MatrixXd>
const llt(lu_m);
6138 if (llt.info() != Eigen::Success) {
6162 double d0 = refpar->parameters()[
Trk::d0];
6163 double z0 = refpar->parameters()[
Trk::z0];
6166 double qoverp = refpar->parameters()[
Trk::qOverP];
6168 if (nperparams > 0) {
6169 d0 += deltaParameters[0];
6170 z0 += deltaParameters[1];
6171 phi += deltaParameters[2];
6172 theta += deltaParameters[3];
6173 qoverp = (trajectory.
m_straightline) ? 0 : .001 * deltaParameters[4] + qoverp;
6185 std::vector < std::pair < double, double >>&scatangles = trajectory.
scatteringAngles();
6186 for (
int i = 0; i < nscat; i++) {
6187 scatangles[i].first += deltaParameters[2 * i + nperparams];
6188 scatangles[i].second += deltaParameters[2 * i + nperparams + 1];
6194 std::vector < double >&delta_ps = trajectory.
brems();
6195 for (
int i = 0; i < nbrem; i++) {
6196 delta_ps[i] += deltaParameters[nperparams + 2 * nscat + i];
6202 std::unique_ptr<const TrackParameters> newper(
6222 const EventContext& evtctx
6233 if (!splitProbContainer.
isValid()) {
6237 std::vector<std::unique_ptr<GXFTrackState>> & states = trajectory.
trackStates();
6244 for (
size_t stateno = 0; stateno < states.size(); stateno++) {
6249 measno += states[stateno-1]->numberOfMeasuredParameters();
6252 std::unique_ptr<GXFTrackState> & state = states[stateno];
6264 const auto *
const rot =
static_cast<const RIO_OnTrack *
>(pMeas);
6265 prd = rot->prepRawData();
6275 const auto &splitProb = splitProbContainer->splitProbability(
pixelCluster);
6276 if (!splitProb.isSplit()) {
6277 ATH_MSG_DEBUG(
"Pixel cluster is not split so no need to update" );
6281 std::unique_ptr < const RIO_OnTrack > newrot;
6282 double *olderror = state->measurementErrors();
6285 double newerror[5] = {-1,-1,-1,-1,-1};
6286 double newres[2] = {-1,-1};
6288 newrot.reset(
m_ROTcreator->correct(*prd, *trackpars, evtctx));
6293 const Amg::MatrixX & covmat = newrot->localCovariance();
6295 newerror[0] = std::sqrt(covmat(0, 0));
6296 newres[0] = newrot->localParameters()[
Trk::locX] - trackpars->parameters()[
Trk::locX];
6297 newerror[1] = std::sqrt(covmat(1, 1));
6298 newres[1] = newrot->localParameters()[
Trk::locY] - trackpars->parameters()[
Trk::locY];
6300 if (
a.cols() != nfitpars) {
6305 for(
int k =0; k<2; k++ ){
6306 const double oldres =
res[measno+k];
6307 res[measno+k] = newres[k];
6308 err[measno+k] = newerror[k];
6310 for (
int i = 0; i < nfitpars; i++) {
6311 if (weightderiv(measno+k, i) == 0) {
6315 b[i] -= weightderiv(measno+k, i) * (oldres / olderror[k] - (newres[k] * olderror[k]) / (newerror[k] * newerror[k]));
6317 for (
int j = i; j < nfitpars; j++) {
6321 weightderiv(measno+k, i) *
6322 weightderiv(measno+k, j) *
6323 ((olderror[k] * olderror[k]) / (newerror[k] * newerror[k]) - 1)
6327 weightderiv(measno+k, i) *= olderror[k] / newerror[k];
6331 state->setMeasurement(std::move(newrot));
6332 state->setMeasurementErrors(newerror);
6347 const EventContext& ctx
6355 std::vector<std::unique_ptr<GXFTrackState>> & states = trajectory.
trackStates();
6361 if (
a.cols() != nfitpars) {
6369 bool outlierremoved =
false;
6370 bool hitrecalibrated =
false;
6372 for (
int stateno = 0; stateno < (int) states.size(); stateno++) {
6373 std::unique_ptr<GXFTrackState> & state = states[stateno];
6381 std::abs(state->trackParameters()->parameters()[
Trk::driftRadius]) > 1.05 * state->associatedSurface().bounds().r()
6386 outlierremoved =
true;
6388 double *errors = state->measurementErrors();
6389 const double olderror = errors[0];
6393 for (
int i = 0; i < nfitpars; i++) {
6394 if (weightderiv(measno, i) == 0) {
6398 b[i] -=
res[measno] * weightderiv(measno, i) / olderror;
6400 for (
int j = i; j < nfitpars; j++) {
6403 a(i, j) - weightderiv(measno, i) * weightderiv(measno, j)
6406 weightderiv(measno, i) = 0;
6410 }
else if (trtrecal) {
6411 double *errors = state->measurementErrors();
6412 const double olderror = errors[0];
6414 const auto *
const thisMeasurement{state->measurement()};
6423 const double trackradius = state->trackParameters()->parameters()[
Trk::driftRadius];
6425 std::unique_ptr<const Trk::RIO_OnTrack> newrot =
nullptr;
6426 const double distance = std::abs(std::abs(trackradius) - dcradius);
6428 if (distance < scalefactor * dcerror && (olderror > 1. || trackradius * oldradius < 0)) {
6430 }
else if (distance > scalefactor * dcerror && olderror < 1.) {
6434 if (newrot !=
nullptr) {
6436 hitrecalibrated =
true;
6440 if ((measno < 0) or (measno >= (
int)
res.size())) {
6441 throw std::runtime_error(
6442 "'res' array index out of range in TrkGlobalChi2Fitter/src/GlobalChi2Fitter.cxx:" + std::to_string(__LINE__)
6446 const double oldres =
res[measno];
6447 const double newres = newradius - state->trackParameters()->parameters()[
Trk::driftRadius];
6448 errors[0] = newerror;
6449 state->setMeasurement(std::move(newrot));
6453 for (
int i = 0; i < nfitpars; i++) {
6454 if (weightderiv(measno, i) == 0) {
6458 b[i] -= weightderiv(measno, i) * (oldres / olderror - (newres * olderror) / (newerror * newerror));
6460 for (
int j = i; j < nfitpars; j++) {
6467 i < nperpars + 2 * nscats &&
6468 (i - nperpars) % 2 == 0
6475 a(i, j) + weightderiv(measno, i) * weightderiv(measno, j) * ((olderror * olderror) / (newerror * newerror) - 1) * weight
6478 weightderiv(measno, i) *= olderror / newerror;
6481 res[measno] = newres;
6482 err[measno] = newerror;
6491 measno += state->numberOfMeasuredParameters();
6495 if (trajectory.
nDOF() < 0) {
6500 if (outlierremoved || hitrecalibrated) {
6509 const EventContext& ctx,
6517 bool trackok =
false;
6519 std::unique_ptr < GXFTrajectory > cleanup_oldtrajectory;
6521 std::unique_ptr < GXFTrajectory > cleanup_newtrajectory;
6527 while (!trackok && oldtrajectory->
nDOF() > 0) {
6529 std::vector<std::unique_ptr<GXFTrackState>> & states = oldtrajectory->
trackStates();
6537 if (nhits != nsihits) {
6541 double maxsipull = -1;
6543 int hitno_maxsipull = -1;
6544 int measno_maxsipull = -1;
6545 int stateno_maxsipull = 0;
6557 for (
int stateno = 0; stateno < (int) states.size(); stateno++) {
6558 std::unique_ptr<GXFTrackState> & state = states[stateno];
6564 double *errors = state->measurementErrors();
6566 const Amg::MatrixX & hitcov = state->measurement()->localCovariance();
6567 const double sinstereo = state->sinStereo();
6568 const double cosstereo = (sinstereo == 0) ? 1 : std::sqrt(1 - sinstereo * sinstereo);
6569 double weight1 = -1;
6571 if (hitcov(0, 0) > trackcov(0, 0)) {
6572 if (sinstereo == 0) {
6573 weight1 = errors[0] * errors[0] - trackcov(0, 0);
6575 weight1 = errors[0] * errors[0] - (
6576 trackcov(0, 0) * cosstereo * cosstereo + 2 *
6577 trackcov(1, 0) * cosstereo * sinstereo + trackcov(1, 1) * sinstereo * sinstereo
6582 const double weight2 = (
6584 errors[1] * errors[1] - trackcov(1, 1) :
6588 double sipull1 = weight1 > 0 ? std::abs(
res[measno] / std::sqrt(weight1)) : -1;
6589 const double sipull2 = (
6591 std::abs(
res[measno + 1] / std::sqrt(weight2)) :
6594 sipull1 = std::max(sipull1, sipull2);
6596 if (sipull1 > maxsipull) {
6597 maxsipull = sipull1;
6598 measno_maxsipull = measno;
6599 state_maxsipull = state.get();
6600 stateno_maxsipull = stateno;
6601 hitno_maxsipull = hitno;
6612 measno += state->numberOfMeasuredParameters();
6616 const double maxpull = maxsipull;
6618 ATH_MSG_DEBUG(
" maxsipull: " << maxsipull <<
" hitno_maxsipull: " <<
6619 hitno_maxsipull <<
" n3sigma: " << n3sigma <<
" cut: " << cut <<
" cut2: " << cut2);
6630 state_maxsipull = oldtrajectory->
trackStates()[stateno_maxsipull].get();
6633 const auto *
const rot =
static_cast<const RIO_OnTrack *
>(pMeas);
6634 prd = rot->prepRawData();
6636 std::unique_ptr < const RIO_OnTrack > broadrot;
6641 Amg::VectorX parameterVector = trackpar_maxsipull->parameters();
6642 const std::unique_ptr<const TrackParameters> trackparForCorrect(
6651 state_maxsipull->trackCovariance())
6655 newerror[0] = newerror[1] = newerror[2] = newerror[3] = newerror[4] = -1;
6656 double newpull = -1;
6657 double newpull1 = -1;
6658 double newpull2 = -1;
6659 double newres1 = -1;
6660 double newres2 = -1;
6661 double newsinstereo = 0;
6674 const Amg::MatrixX & covmat = broadrot->localCovariance();
6676 if (state_maxsipull->
sinStereo() != 0) {
6677 const auto [covEigenValueSmall, covStereoAngle] = principalComponentAnalysis2x2(covmat);
6678 newerror[0] = std::sqrt(covEigenValueSmall);
6679 newsinstereo = std::sin(covStereoAngle);
6681 newerror[0] = std::sqrt(covmat(0, 0));
6684 const double cosstereo = (newsinstereo == 0) ? 1. : std::sqrt(1 - newsinstereo * newsinstereo);
6686 if (cosstereo != 1.) {
6688 cosstereo * (broadrot->localParameters()[
Trk::locX] - trackpar_maxsipull->parameters()[
Trk::locX]) +
6689 newsinstereo * (broadrot->localParameters()[
Trk::locY] - trackpar_maxsipull->parameters()[
Trk::locY])
6692 newres1 = broadrot->localParameters()[
Trk::locX] - trackpar_maxsipull->parameters()[
Trk::locX];
6695 if (newerror[0] == 0.0) {
6696 ATH_MSG_WARNING(
"Measurement error is zero or negative, treating as outlier");
6699 newpull1 = std::abs(newres1 / newerror[0]);
6703 newerror[1] = std::sqrt(covmat(1, 1));
6704 newres2 = broadrot->localParameters()[
Trk::locY] - trackpar_maxsipull->parameters()[
Trk::locY];
6705 newpull2 = std::abs(newres2 / newerror[1]);
6708 newpull = std::max(newpull1, newpull2);
6714 (newerror[0] > 1.5 * olderror[0] || newerror[1] > 1.5 * std::abs(olderror[1]))
6716 if ((measno_maxsipull < 0) or(measno_maxsipull >= (
int)
res.size())) {
6717 throw std::runtime_error(
6718 "'res' array index out of range in TrkGlobalChi2Fitter/src/GlobalChi2Fitter.cxx:" + std::to_string(__LINE__)
6723 newtrajectory = oldtrajectory;
6725 if (
a.cols() != nfitpars) {
6729 const double oldres1 =
res[measno_maxsipull];
6730 res[measno_maxsipull] = newres1;
6731 err[measno_maxsipull] = newerror[0];
6733 for (
int i = 0; i < nfitpars; i++) {
6734 if (weightderiv(measno_maxsipull, i) == 0) {
6738 b[i] -= weightderiv(measno_maxsipull, i) * (oldres1 / olderror[0] - (newres1 * olderror[0]) / (newerror[0] * newerror[0]));
6740 for (
int j = i; j < nfitpars; j++) {
6744 weightderiv(measno_maxsipull, i) *
6745 weightderiv(measno_maxsipull, j) *
6746 ((olderror[0] * olderror[0]) / (newerror[0] * newerror[0]) - 1)
6750 weightderiv(measno_maxsipull, i) *= olderror[0] / newerror[0];
6754 const double oldres2 =
res[measno_maxsipull + 1];
6755 res[measno_maxsipull + 1] = newres2;
6756 err[measno_maxsipull + 1] = newerror[1];
6758 for (
int i = 0; i < nfitpars; i++) {
6759 if (weightderiv(measno_maxsipull + 1, i) == 0) {
6763 b[i] -= weightderiv(measno_maxsipull + 1, i) * (oldres2 / olderror[1] - (newres2 * olderror[1]) / (newerror[1] * newerror[1]));
6765 for (
int j = i; j < nfitpars; j++) {
6769 weightderiv(measno_maxsipull + 1, i) *
6770 weightderiv(measno_maxsipull + 1, j) *
6771 ((olderror[1] * olderror[1]) / (newerror[1] * newerror[1]) - 1)
6776 weightderiv(measno_maxsipull + 1, i) *= olderror[1] / newerror[1];
6781 "Recovering outlier, hitno=" << hitno_maxsipull <<
" measno=" <<
6782 measno_maxsipull <<
" pull=" << maxsipull <<
" olderror_0=" <<
6783 olderror[0] <<
" newerror_0=" << newerror[0] <<
" olderror_1=" <<
6784 olderror[1] <<
" newerror_1=" << newerror[1]
6793 ((n3sigma < 2 && maxsipull > cut2 && maxsipull < cut) || n3sigma > 1) &&
6803 "Removing outlier, hitno=" << hitno_maxsipull <<
", measno=" <<
6804 measno_maxsipull <<
" pull=" << maxsipull
6811 cleanup_newtrajectory = std::make_unique<GXFTrajectory>(*oldtrajectory);
6812 newtrajectory = cleanup_newtrajectory.get();
6814 if (newa.cols() != nfitpars) {
6820 if ((measno_maxsipull < 0) or(measno_maxsipull >= (
int)
res.size())) {
6821 throw std::runtime_error(
6822 "'res' array index out of range in TrkGlobalChi2Fitter/src/GlobalChi2Fitter.cxx:" + std::to_string(__LINE__)
6826 const double oldres1 =
res[measno_maxsipull];
6827 newres[measno_maxsipull] = 0;
6829 for (
int i = 0; i < nfitpars; i++) {
6830 if (weightderiv(measno_maxsipull, i) == 0) {
6834 newb[i] -= weightderiv(measno_maxsipull, i) * oldres1 / olderror[0];
6836 for (
int j = i; j < nfitpars; j++) {
6840 weightderiv(measno_maxsipull, i) *
6841 weightderiv(measno_maxsipull, j)
6845 newweightderiv(measno_maxsipull, i) = 0;
6849 const double oldres2 =
res[measno_maxsipull + 1];
6850 newres[measno_maxsipull + 1] = 0;
6852 for (
int i = 0; i < nfitpars; i++) {
6853 if (weightderiv(measno_maxsipull + 1, i) == 0) {
6857 newb[i] -= weightderiv(measno_maxsipull + 1, i) * oldres2 / olderror[1];
6859 for (
int j = i; j < nfitpars; j++) {
6860 if (weightderiv(measno_maxsipull + 1, j) == 0) {
6867 weightderiv(measno_maxsipull + 1, i) *
6868 weightderiv(measno_maxsipull + 1, j)
6872 newweightderiv(measno_maxsipull + 1, i) = 0;
6876 newtrajectory->
setOutlier(stateno_maxsipull);
6890 for (
int it = 0; it <
m_maxit; ++it) {
6900 ctx, cache, *newtrajectory, it, *newap, *newbp, lu_m, doderiv);
6916 const double oldchi2 = oldtrajectory->
chi2() / oldtrajectory->
nDOF();
6917 const double newchi2 = (newtrajectory->
nDOF() > 0) ? newtrajectory->
chi2() / newtrajectory->
nDOF() : 0;
6920 if (newtrajectory->
nDOF() != oldtrajectory->
nDOF() && maxsipull > cut2) {
6921 mindiff = (oldchi2 > .33 *
m_chi2cut || noutl > 0) ? .8 : 1.;
6923 if (noutl == 0 && maxsipull < cut - .5 && oldchi2 < .5 *
m_chi2cut) {
6928 if (newchi2 > oldchi2 || (newchi2 > oldchi2 - mindiff && newchi2 > .33 * oldchi2)) {
6929 ATH_MSG_DEBUG(
"Outlier not confirmed, keeping old trajectory");
6937 (void)cleanup_oldtrajectory.release();
6938 return oldtrajectory;
6940 if (oldtrajectory != newtrajectory) {
6941 cleanup_oldtrajectory = std::move(cleanup_newtrajectory);
6942 oldtrajectory = newtrajectory;
6949 Eigen::LLT < Eigen::MatrixXd >
const lltOfW(
a);
6950 if (lltOfW.info() == Eigen::Success) {
6954 const int ncols =
a.cols();
6955 Amg::MatrixX const weightInvAMG = Amg::MatrixX::Identity(ncols, ncols);
6956 fullcov = lltOfW.solve(weightInvAMG);
6974 oldtrajectory->
nDOF() > 0 &&
6983 (void)cleanup_oldtrajectory.release();
6984 return oldtrajectory;
6996 if (
const auto *pMeas{hit->measurement()};
7002 nrealmeas += hit->numberOfMeasuredParameters();
7012 if (
const auto *pMeas{hit->measurement()};
7018 for (
int i = measindex; i < measindex + hit->numberOfMeasuredParameters(); i++) {
7020 cache.
m_derivmat(i, j) = derivs(measindex2, j) * errors[measindex2];
7021 if ((j == 4 && !oldtrajectory.
m_straightline) || j >= nperpars + 2 * nscat) {
7029 measindex += hit->numberOfMeasuredParameters();
7030 }
else if (hit->materialEffects() ==
nullptr) {
7031 measindex2 += hit->numberOfMeasuredParameters();
7037 const EventContext & ctx,
7045 std::unique_ptr<const TrackParameters> per(
nullptr);
7048 std::unique_ptr<const TrackParameters> prevpar(
7053 std::vector<std::pair<const Layer *, const Layer *>> & upstreamlayers = oldtrajectory.
upstreamMaterialLayers();
7056 for (
const auto & [layer1, layer2] : upstreamlayers | std::views::reverse) {
7057 if (prevpar ==
nullptr) {
7062 const Layer *layer = layer1 !=
nullptr ? layer1 : layer2;
7064 const DistanceSolution distsol = layer->surfaceRepresentation().straightLineDistanceEstimate(
7065 prevpar->position(), prevpar->momentum().unit()
7067 const double distance = getDistance(distsol);
7070 if (std::abs(distance) < 0.01) {
7074 if (distsol.
first() * distsol.
second() < 0 && !first) {
7079 if (first && distance > 0) {
7083 std::unique_ptr<const TrackParameters> layerpar(
7087 layer->surfaceRepresentation(),
7095 if (layerpar ==
nullptr) {
7099 if (layer->surfaceRepresentation().bounds().inside(layerpar->localPosition())) {
7103 prevpar = std::move(layerpar);
7116 if (startfactor > 0.5) {
7117 std::unique_ptr<const TrackParameters> updatedpar =
m_matupdator->update(
7121 if (updatedpar !=
nullptr) {
7141 if (endfactor > 0.5) {
7142 std::unique_ptr<const TrackParameters> updatedpar =
m_matupdator->update(
7146 if (updatedpar !=
nullptr) {
7152 if (prevpar !=
nullptr) {
7164 if (per ==
nullptr) {
7165 ATH_MSG_DEBUG(
"Failed to extrapolate to perigee, returning 0");
7184 std::unique_ptr<GXFTrackState>
7186 const EventContext & ctx,
7193 if (per ==
nullptr) {
7197 ATH_MSG_DEBUG(
"Final perigee: " << *per <<
" pos: " << per->position() <<
" pT: " << per->pT());
7203 const std::vector<std::unique_ptr<TrackParameters>> & hc,
7204 std::set<Identifier> & id_set,
7205 std::set<Identifier> & sct_set,
7215 for (
const std::unique_ptr<TrackParameters> & tp : hc) {
7221 if (tp ==
nullptr) {
7230 const TrkDetElementBase * de = tp->associatedSurface().associatedDetectorElement();
7232 if (de ==
nullptr) {
7243 if (id_set.find(
id) != id_set.end()) {
7262 }
else if (
m_DetID->is_sct(
id)) {
7273 }
else if (
m_DetID->is_sct(
id)) {
7282 const Identifier os = e->otherSide()->identify();
7292 if (sct_set.find(os) != sct_set.end()) {
7293 ++rv.m_sct_double_hole;
7322 for (
const std::unique_ptr<GXFTrackState> & s : trajectory.
trackStates()) {
7334 std::vector<std::reference_wrapper<GXFTrackState>> rv;
7342 for (
const std::unique_ptr<GXFTrackState> & s : trajectory.
trackStates()) {
7356 rv.emplace_back(*s);
7363 const TrkDetElementBase * de = s->trackParameters()->associatedSurface().associatedDetectorElement();
7365 if (de !=
nullptr) {
7378 if (s.get() == lastmeas) {
7388 const EventContext & ctx,
7389 const std::vector<std::reference_wrapper<GXFTrackState>> & states
7401 constexpr uint min_meas = 3;
7402 if (std::count_if(states.begin(), states.end(), [](
const GXFTrackState & s){ return s.getStateType(TrackStateOnSurface::Measurement); }) < min_meas) {
7406 bool seen_meas =
false;
7408 std::set<Identifier> id_set;
7409 std::set<Identifier> sct_set;
7415 for (std::size_t i = 0; i < states.size() - 1; i++) {
7438 const double dist = (beg.trackParameters()->position() - end.trackParameters()->position()).norm();
7440 const bool zStartValid = std::abs(beg.trackParameters()->position().z())<10000.;
7442 ATH_MSG_DEBUG(
"Pathological track parameter well outside of detector");
7443 ATH_MSG_DEBUG(
"Propagator might have issue with this, skipping");
7444 ATH_MSG_VERBOSE(
"dumping track parameters " << *(beg.trackParameters()));
7453 if (seen_meas && dist >= 2.5 && zStartValid) {
7460 std::optional<std::vector<std::unique_ptr<TrackParameters>>> & hc = beg.getHoles();
7461 std::vector<std::unique_ptr<TrackParameters>> states;
7469 if (hc.has_value()) {
7470 states = std::move(*hc);
7508 std::vector<std::unique_ptr<Trk::TrackParameters>>
const bl =
m_extrapolator->extrapolateBlindly(
7529 const EventContext & ctx,
7535 auto trajectory = std::make_unique<Trk::TrackStates>();
7543 std::unique_ptr<GXFTrackState> perigee_ts =
makeTrackFindPerigee(ctx, cache, oldtrajectory, matEffects);
7545 if (perigee_ts ==
nullptr) {
7551 trajectory->reserve(tmptrajectory.
trackStates().size());
7557 hit->resetTrackCovariance();
7563 hit->materialEffects()))
7569 auto trackState = hit->trackStateOnSurface();
7570 hit->resetTrackCovariance();
7571 trajectory->emplace_back(trackState.release());
7574 auto qual = std::make_unique<FitQuality>(tmptrajectory.
chi2(), tmptrajectory.
nDOF());
7594 std::unique_ptr<Track> rv = std::make_unique<Track>(info, std::move(trajectory), std::move(qual));
7604 std::unique_ptr<TrackSummary>
ts = std::make_unique<TrackSummary>();
7613 std::optional<TrackHoleCount> hole_count;
7620 std::vector<std::reference_wrapper<GXFTrackState>>
const states =
holeSearchStates(tmptrajectory);
7638 if (hole_count.has_value()) {
7651 rv->setTrackSummary(std::move(
ts));
7660 const EventContext & ctx,
7669 std::vector<std::unique_ptr<TrackParameters>> rv =
m_extrapolator->extrapolateStepwise(
7683 &rv.front()->associatedSurface() == &src.associatedSurface() ||
7684 trackParametersClose(*rv.front(), src, 0.001) ||
7688 rv.front().reset(
nullptr);
7699 &rv.back()->associatedSurface() == &src.associatedSurface() ||
7700 trackParametersClose(*rv.back(), src, 0.001) ||
7704 rv.back().reset(
nullptr);
7711 const EventContext & ctx,
7719 std::unique_ptr<const TrackParameters> rv;
7720 std::optional<TransportJacobian> jac{};
7731 if (rv !=
nullptr && calcderiv) {
7736 std::optional<std::vector<std::unique_ptr<TrackParameters>>> extrapolation;
7745 std::move(extrapolation)
7750 const EventContext & ctx,
7761 ctx, prev,
ts, propdir, bf, calcderiv, holesearch
7764 if (rv.m_parameters ==
nullptr) {
7765 propdir = invertPropdir(propdir);
7768 ctx, prev,
ts, propdir, bf, calcderiv, holesearch
7776 const EventContext& ctx,
7783 std::vector<std::unique_ptr<GXFTrackState>> & states = trajectory.
trackStates();
7786 std::unique_ptr<const TrackParameters> tmptrackpar;
7788 for (
int hitno = nstatesupstream - 1; hitno >= 0; hitno--) {
7789 const Surface &surf1 = states[hitno]->associatedSurface();
7796 const double distance = getDistance(distsol);
7818 (rv.m_parameters !=
nullptr) &&
7819 (prevtrackpar->
position() - rv.m_parameters->position()).mag() > 5 * mm
7825 if (rv.m_parameters ==
nullptr) {
7826 ATH_MSG_DEBUG(
"propagation failed, prev par: " << *prevtrackpar <<
7827 " pos: " << prevtrackpar->
position() <<
" destination surface: " << surf1);
7831 states[hitno]->setTrackParameters(std::move(rv.m_parameters));
7832 const TrackParameters *currenttrackpar = states[hitno]->trackParameters();
7833 const Surface &surf = states[hitno]->associatedSurface();
7835 if (rv.m_jacobian != std::nullopt) {
7837 states[hitno]->materialEffects() !=
nullptr &&
7838 states[hitno]->materialEffects()->deltaE() != 0 &&
7839 states[hitno]->materialEffects()->sigmaDeltaE() <= 0 &&
7842 const double p = 1. / std::abs(currenttrackpar->parameters()[
Trk::qOverP]);
7843 const double de = std::abs(states[hitno]->materialEffects()->deltaE());
7844 const double mass = trajectory.
mass();
7845 const double newp = std::sqrt(p * p + 2 * de * std::sqrt(mass * mass + p * p) + de * de);
7846 (*rv.m_jacobian) (4, 4) = ((p + p * de / std::sqrt(p * p + mass * mass)) / newp) * p * p / (newp * newp);
7849 states[hitno]->setJacobian(*rv.m_jacobian);
7850 }
else if (calcderiv) {
7857 if (meff !=
nullptr && hitno != 0) {
7859 surf, *meff, *states[hitno]->trackParameters(), trajectory.
mass(), -1
7862 if (std::holds_alternative<FitterStatusCode>(
r)) {
7863 return std::get<FitterStatusCode>(
r);
7866 tmptrackpar = std::move(std::get<std::unique_ptr<const TrackParameters>>(
r));
7867 prevtrackpar = tmptrackpar.get();
7869 prevtrackpar = currenttrackpar;
7875 for (
int hitno = nstatesupstream; hitno < (int) states.size(); hitno++) {
7876 const Surface &surf = states[hitno]->associatedSurface();
7880 const double distance = getDistance(distsol);
7897 (rv.m_parameters !=
nullptr) &&
7899 (prevtrackpar->
position() - rv.m_parameters->position()).mag() > 5 * mm
7904 if (rv.m_parameters ==
nullptr) {
7905 ATH_MSG_DEBUG(
"propagation failed, prev par: " << *prevtrackpar <<
7906 " pos: " << prevtrackpar->
7907 position() <<
" destination surface: " << surf);
7911 if (rv.m_jacobian != std::nullopt) {
7913 states[hitno]->materialEffects() !=
nullptr &&
7914 states[hitno]->materialEffects()->deltaE() != 0 &&
7915 states[hitno]->materialEffects()->sigmaDeltaE() <= 0 &&
7918 const double p = 1 / std::abs(rv.m_parameters->parameters()[
Trk::qOverP]);
7919 const double de = std::abs(states[hitno]->materialEffects()->deltaE());
7920 const double mass = trajectory.
mass();
7921 double newp = p * p - 2 * de * std::sqrt(mass * mass + p * p) + de * de;
7924 newp = std::sqrt(newp);
7927 (*rv.m_jacobian) (4, 4) = ((p - p * de / std::sqrt(p * p + mass * mass)) / newp) * p * p / (newp * newp);
7930 states[hitno]->setJacobian(*rv.m_jacobian);
7931 }
else if (calcderiv) {
7938 if (meff !=
nullptr) {
7940 surf, *meff, *rv.m_parameters, trajectory.
mass(), +1
7943 if (std::holds_alternative<FitterStatusCode>(
r)) {
7944 return std::get<FitterStatusCode>(
r);
7947 rv.m_parameters = std::move(std::get<std::unique_ptr<const TrackParameters>>(
r));
7950 states[hitno]->setTrackParameters(std::move(rv.m_parameters));
7951 prevtrackpar = states[hitno]->trackParameters();
7964 const AmgVector(5) & old = param.parameters();
7970 ATH_MSG_DEBUG(
"Angles out of range, phi: " << newphi <<
" theta: " << newtheta);
7974 double newqoverp = 0;
7980 const double oldp = std::abs(1 / old[
Trk::qOverP]);
7981 const double newp2 = oldp * oldp -
sign * 2 * std::abs(meff.
deltaE()) * std::sqrt(mass * mass + oldp * oldp) + meff.
deltaE() * meff.
deltaE();
7988 newqoverp = std::copysign(1 / std::sqrt(newp2), old[
Trk::qOverP]);
7995 old[0], old[1], newphi, newtheta, newqoverp, std::nullopt
8007 using Matrix55 = Eigen::Matrix<double, 5, 5>;
8009 Matrix55 initialjac;
8010 initialjac.setZero();
8011 initialjac(4, 4) = 1;
8013 Matrix55 jacvertex(initialjac);
8015 std::vector<Matrix55, Eigen::aligned_allocator<Matrix55>> jacscat(trajectory.
numberOfScatterers(), initialjac);
8016 std::vector<Matrix55, Eigen::aligned_allocator<Matrix55>> jacbrem(trajectory.
numberOfBrems(), initialjac);
8018 std::vector<std::unique_ptr<GXFTrackState>> & states = trajectory.
trackStates();
8027 for (
const bool forward : {
false,
true}) {
8029 hit_begin = nstatesupstream;
8030 hit_end = (int) states.size();
8031 scatno = nscatupstream;
8032 bremno = nbremupstream;
8034 hit_begin = nstatesupstream - 1;
8041 int hitno = hit_begin;
8042 forward ? (hitno < hit_end) : (hitno >= hit_end);
8043 hitno += (forward ? 1 : -1)
8046 state = states[hitno].get();
8050 if (fillderivmat && state->
derivatives().cols() != nfitpars) {
8058 const int jmaxbrem = 4;
8060 if (hitno == (forward ? hit_end - 1 : 0)) {
8061 if (!fillderivmat) {
8069 Eigen::Matrix<double, 5, 5> & jac = state->
jacobian();
8071 if (hitno == nstatesupstream + (forward ? 0 : -1)) {
8072 jacvertex.block<4, 5>(0, 0) = jac.block<4, 5>(0, 0);
8073 jacvertex(4, 4) = jac(4, 4);
8083 jcnt = jmax - jmin + 1;
8085 lp_bgn = forward ? nscatupstream : nscatupstream - 1;
8088 for (
int i = lp_bgn; forward ? (i < lp_end) : (i > lp_end); i += (forward ? 1 : -1)) {
8090 i == scatno + (forward ? -1 : 1) &&
8091 prevstate !=
nullptr &&
8095 jacscat[i].block(0, jmin, 4, jcnt) = jac.block(0, jmin, 4, jcnt);
8096 jacscat[i](4, 4) = jac(4, 4);
8098 calculateJac(jac, jacscat[i], jmin, jmax);
8102 Eigen::MatrixXd & derivmat = state->
derivatives();
8103 const int scatterPos = nperpars + 2 * i;
8105 derivmat.block<4, 2>(0, scatterPos) = (forward ? 1 : -1) * jacscat[i].block<4, 2>(0, 2);
8111 jcnt = jmax - jmin + 1;
8113 lp_bgn = forward ? nbremupstream : nbremupstream - 1;
8116 for (
int i = lp_bgn; forward ? (i < lp_end) : (i > lp_end); i += (forward ? 1 : -1)) {
8118 i == bremno + (forward ? -1 : 1) &&
8123 jacbrem[i].block(0, jmin, 4, jcnt) = jac.block(0, jmin, 4, jcnt);
8124 jacbrem[i](4, 4) = jac(4, 4);
8126 calculateJac(jac, jacbrem[i], jmin, jmax);
8130 Eigen::MatrixXd & derivmat = state->
derivatives();
8131 const int scatterPos = nperpars + 2 * nscats + i;
8133 derivmat.block<5, 1>(0, scatterPos) = (forward ? .001 : -.001) * jacbrem[i].block<5, 1>(0, 4);
8137 calculateJac(jac, jacvertex, 0, 4);
8141 Eigen::MatrixXd & derivmat = state->
derivatives();
8142 derivmat.block(0, 0, 4, nperpars) = jacvertex.block(0, 0, 4, nperpars);
8144 if (nperpars == 5) {
8145 derivmat.col(4).segment(0, 4) *= .001;
8146 derivmat(4, 4) = .001 * jacvertex(4, 4);
8152 (!trajectory.
prefit() || states[hitno]->materialEffects()->deltaE() == 0)
8154 scatno += (forward ? 1 : -1);
8158 states[hitno]->materialEffects() &&
8159 states[hitno]->materialEffects()->sigmaDeltaE() > 0
8161 bremno += (forward ? 1 : -1);
8164 prevstate = states[hitno].get();
8173 bool onlylocal)
const {
8179 std::vector<std::unique_ptr<GXFTrackState>> & states = trajectory.
trackStates();
8181 std::vector < int >
indices(states.size());
8183 int i = nstatesupstream;
8184 for (
int j = 0; j < (int) states.size(); j++) {
8185 if (j < nstatesupstream) {
8192 for (
int stateno = 0; stateno < (int) states.size(); stateno++) {
8193 if (stateno == 0 || stateno == nstatesupstream) {
8194 prevstate =
nullptr;
8197 std::unique_ptr<GXFTrackState> & state = states[
index];
8198 if (state->materialEffects() !=
nullptr) {
8199 prevstate = state.get();
8203 if (!state->hasTrackCovariance()) {
8204 state->zeroTrackCovariance();
8206 AmgMatrix(5, 5) & trackerrmat = state->trackCovariance();
8208 if ((prevstate !=
nullptr) &&
8212 Eigen::Matrix<double, 5, 5> & jac = state->jacobian();
8213 const AmgMatrix(5, 5)& prevcov = states[
indices[stateno - 1]]->trackCovariance();
8215 trackerrmat = jac * prevcov * jac.transpose();
8219 trackerrmat = derivatives * fullcovmat * derivatives.transpose();
8229 bool errorok =
true;
8230 for (
int i = 0; i < 5; i++) {
8233 && trackerrmat(i, i) > meascov(j, j)) {
8235 const double scale = std::sqrt(meascov(j, j) / trackerrmat(i, i));
8236 trackerrmat(i, i) = meascov(j, j);
8237 for (
int k = 0; k < 5; k++) {
8239 trackerrmat(k, i) *= scale;
8247 for (
int i = 0; i < 5; i++) {
8251 for (
int j = 0; j < 5; j++) {
8259 trackerrmat(4, 4) = 1e-20;
8263 state->trackParameters();
8265 std::optional<
AmgMatrix(5, 5)> trkerrmat;
8267 if (state->hasTrackCovariance()) {
8268 trkerrmat = (state->trackCovariance());
8270 trkerrmat = std::nullopt;
8273 const AmgVector(5) & tpars = tmptrackpar->parameters();
8274 std::unique_ptr<const TrackParameters> trackpar(
8280 std::move(trkerrmat))
8282 state->setTrackParameters(std::move(trackpar));
8285 if (errorok && trajectory.
nDOF() > 0) {
8286 fitQual =
m_updator->fullStateFitQuality(
8287 *state->trackParameters(),
8295 state->setFitQuality(fitQual);
8297 prevstate = state.get();
8301 std::optional<TransportJacobian>
8303 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) {
8339 if (i == 0 && thiscylsurf) {
8340 vecminuseps[i] = -std::remainder(-vecminuseps[i], 2 *
M_PI * previousSurface.
bounds().
r());
8341 }
else if (i == 1 && thisdiscsurf) {
8342 vecpluseps[i] = -std::remainder(-vecpluseps[i], 2 *
M_PI);
8347 std::unique_ptr<const TrackParameters> parpluseps(
8357 const 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++) {
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){
8480 testrot = &testcrot->rioOnTrack(0);
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){
8493 testrot = &testcrot->rioOnTrack(0);
8500 (testrot !=
nullptr) &&
8517 if (cond_obj ==
nullptr) {
8526 std::stringstream
msg;
8528 throw std::runtime_error(
msg.str());
8535 if (geometry !=
nullptr) {
8536 cache.
m_caloEntrance = geometry->trackingVolume(
"InDet::Containers::InnerDetector");
8556 if (geometry !=
nullptr) {
8557 cache.
m_msEntrance = geometry->trackingVolume(
"MuonSpectrometerEntrance");
Scalar perp() const
perp method - perpendicular length
Scalar deltaPhi(const MatrixBase< Derived > &vec) const
Scalar mag() const
mag method
#define ATH_CHECK
Evaluate an expression and check for errors.
#define ATH_MSG_VERBOSE(x)
#define ATH_MSG_WARNING(x)
This class provides an interface to generate or decode an identifier for the upper levels of the dete...
double charge(const T &p)
std::vector< size_t > vec
#define AmgSymMatrix(dim)
#define AmgMatrix(rows, cols)
std::pair< std::vector< unsigned int >, bool > res
static const uint32_t nHits
void diff(const Jet &rJet1, const Jet &rJet2, std::map< std::string, double > varDiff)
Difference between jets - Non-Class function required by trigger.
void getInitializedCache(MagField::AtlasFieldCache &cache) const
get B field cache for evaluation as a function of 2-d or 3-d position.
DataModel_detail::const_iterator< DataVector > const_iterator
Standard const_iterator.
const_iterator end() const noexcept
Return a const_iterator pointing past the end of the collection.
const_iterator begin() const noexcept
Return a const_iterator pointing at the beginning of the collection.
size_type size() const noexcept
Returns the number of elements in the collection.
bool empty() const noexcept
Returns true if the collection is empty.
bool is_valid() const
Check if id is in a valid state.
Class to hold geometrical description of a silicon detector element.
bool solenoidOn() const
status of the magnets
void getFieldZR(const double *ATH_RESTRICT xyz, double *ATH_RESTRICT bxyz, double *ATH_RESTRICT deriv=nullptr)
get B field valaue on the z-r plane at given position works only inside the solenoid.
virtual bool isValid() override final
Can the handle be successfully dereferenced?
virtual std::span< T *const > arrayObjects()=0
Return all objects of the Array non-const we can still modify the T.
Base class for all CompetingRIOsOnTack implementations, extends the common MeasurementBase.
virtual const RIO_OnTrack & rioOnTrack(unsigned int) const =0
returns the RIO_OnTrack (also known as ROT) objects depending on the integer.
Bounds for a cylindrical Surface.
virtual double r() const override final
This method returns the radius.
double halflengthZ() const
This method returns the halflengthZ.
Class to describe a cylindrical detector layer for tracking, it inhertis from both,...
Class for a CylinderSurface in the ATLAS detector.
virtual const CylinderBounds & bounds() const override final
This method returns the CylinderBounds by reference (NoBounds is not possible for cylinder)
Bounds for a cylindrical Volume, the decomposeToSurfaces method creates a vector of up to 6 surfaces:
Class to describe the bounds for a planar DiscSurface.
double rMax() const
This method returns outer radius.
double rMin() const
This method returns inner radius.
Class to describe a disc-like detector layer for tracking, it inhertis from both, Layer base class an...
Class for a DiscSurface in the ATLAS detector.
const SurfaceBounds & bounds() const override final
This method returns the bounds by reference.
Access to distance solutions.
double second() const
Distance to second intersection solution along direction (for a cylinder surface)
int numberOfSolutions() const
Number of intersection solutions.
double first() const
Distance to first intersection solution along direction.
double sigmaDeltaE() const
returns the symmatric error
double deltaE() const
returns the
int numberDoF() const
returns the number of degrees of freedom of the overall track or vertex fit as integer
double chiSquared() const
returns the of the overall track fit
Status codes for track fitters.
@ ExtrapolationFailureDueToSmallMomentum
extrapolation failed due to small momentum
@ OutlierLogicFailure
outlier logic failed
@ ExtrapolationFailure
extrapolation failed
class that is similar to MaterialEffectsOnTrack, but has 'set' methods for more flexibility during tr...
void setSigmaDeltaE(double)
double deltaTheta() const
double sigmaDeltaEPos() const
double sigmaDeltaE() const
const Surface & associatedSurface() const
double sigmaDeltaPhi() const
void setScatteringSigmas(double, double)
double sigmaDeltaTheta() const
void setMaterialProperties(const MaterialProperties *)
Set the material properties of this material effects instance.
const MaterialProperties * materialProperties() const
bool getStateType(TrackStateOnSurface::TrackStateOnSurfaceType type) const
Retrieve the value of a specific type bit.
Amg::MatrixX & derivatives()
Eigen::Matrix< double, 5, 5 > & jacobian()
void setSinStereo(double)
const MeasurementBase * measurement(void)
TrackState::MeasurementType measurementType()
bool hasTrackCovariance(void) const
void setMeasurement(std::unique_ptr< const MeasurementBase >)
void setRecalibrated(bool)
void setMeasurementErrors(const double *)
const TrackParameters * trackParameters(void) const
double * measurementErrors()
GXFMaterialEffects * materialEffects()
void setTrackParameters(std::unique_ptr< const TrackParameters >)
const Surface & associatedSurface() const
Internal representation of the track, used in the track fit.
int numberOfOutliers() const
std::vector< double > & brems()
void setNumberOfScatterers(int)
int numberOfSiliconHits() const
std::vector< std::pair< const Layer *, const Layer * > > & upstreamMaterialLayers()
int numberOfFitParameters() const
int numberOfTRTPrecHits() const
Amg::MatrixX & weightedResidualDerivatives()
const std::vector< std::unique_ptr< GXFTrackState > > & trackStates() const
int numberOfTRTTubeHits() const
int numberOfTRTHits() const
int numberOfUpstreamBrems() const
double totalEnergyLoss() const
int numberOfUpstreamScatterers() const
int numberOfUpstreamStates() const
MagneticFieldProperties m_fieldprop
void setNumberOfPerigeeParameters(int)
const TrackParameters * referenceParameters()
void setNumberOfBrems(int)
std::vector< std::pair< double, double > > & scatteringSigmas()
std::pair< GXFTrackState *, GXFTrackState * > findFirstLastMeasurement(void)
void addMaterialState(std::unique_ptr< GXFTrackState >, int index=-1)
int numberOfScatterers() const
bool addMeasurementState(std::unique_ptr< GXFTrackState >, int index=-1)
std::vector< std::pair< double, double > > & scatteringAngles()
void addBasicState(std::unique_ptr< GXFTrackState >, int index=-1)
Amg::VectorX & residuals()
void setBrems(std::vector< double > &)
int numberOfPerigeeParameters() const
int numberOfBrems() const
void setOutlier(int, bool isoutlier=true)
void updateTRTHitCount(int index, float oldError)
void setScatteringAngles(std::vector< std::pair< double, double > > &)
void setReferenceParameters(std::unique_ptr< const TrackParameters >)
Gaudi::Property< int > m_maxoutliers
ToolHandle< IBoundaryCheckTool > m_boundaryCheckTool
std::unique_ptr< Track > makeTrack(const EventContext &ctx, Cache &, GXFTrajectory &, const ParticleHypothesis) const
void tryToConverge(const Cache &cache, GXFTrajectory &trajectory, const int it) const
void makeProtoState(Cache &, GXFTrajectory &, const TrackStateOnSurface *, int index=-1) const
Gaudi::Property< double > m_p
FitterStatusCode calculateTrackParameters(const EventContext &ctx, GXFTrajectory &, bool) const
Track * myfit(const EventContext &ctx, Cache &, GXFTrajectory &, const TrackParameters &, const RunOutlierRemoval runOutlier=false, const ParticleHypothesis matEffects=nonInteracting) const
void fillResidualsAndErrors(const EventContext &ctx, const Cache &cache, GXFTrajectory &trajectory, const int it, Amg::VectorX &b, int &bremno_maxbrempull, GXFTrackState *&state_maxbrempull) const
Gaudi::Property< double > m_scalefactor
Gaudi::Property< bool > m_rejectLargeNScat
std::variant< std::unique_ptr< const TrackParameters >, FitterStatusCode > updateEnergyLoss(const Surface &, const GXFMaterialEffects &, const TrackParameters &, double, int) const
ToolHandle< Trk::ITrkMaterialProviderTool > m_caloMaterialProvider
Gaudi::Property< bool > m_trtrecal
bool processTrkVolume(Cache &, const Trk::TrackingVolume *tvol) const
ToolHandle< IPropagator > m_propagator
Gaudi::Property< bool > m_straightlineprop
Gaudi::Property< bool > m_fiteloss
std::optional< GlobalChi2Fitter::TrackHoleCount > holeSearchProcess(const EventContext &ctx, const std::vector< std::reference_wrapper< GXFTrackState > > &states) const
Conduct a hole search between a list of states, possibly reusing existing information.
static std::optional< std::pair< Amg::Vector3D, double > > addMaterialFindIntersectionDisc(Cache &cache, const DiscSurface &surface, const TrackParameters ¶m1, const TrackParameters ¶m2, const ParticleHypothesis mat)
Find the intersection of a set of track parameters onto a disc surface.
ToolHandle< INavigator > m_navigator
void updateSystemWithMaxBremPull(GXFTrajectory &trajectory, const int bremno_maxbrempull, GXFTrackState *state_maxbrempull, Amg::SymMatrixX &a) const
SG::ReadHandleKey< Trk::ClusterSplitProbabilityContainer > m_clusterSplitProbContainer
virtual Track * alignmentFit(AlignmentCache &, const Track &, const RunOutlierRemoval runOutlier=false, const ParticleHypothesis matEffects=Trk::nonInteracting) const override
virtual ~GlobalChi2Fitter()
ToolHandle< IRIO_OnTrackCreator > m_broadROTcreator
std::unique_ptr< GXFTrackState > makeTrackFindPerigee(const EventContext &, Cache &, GXFTrajectory &, const ParticleHypothesis) const
FitterStatusCode updateFitParameters(GXFTrajectory &, const Amg::VectorX &, const Amg::SymMatrixX &) const
Method to update peregee parameters, scattering angles, and brems.
void addMaterial(const EventContext &ctx, Cache &, GXFTrajectory &, const TrackParameters *, ParticleHypothesis) const
static void fillFirstLastMeasurement(Cache &cache, GXFTrajectory &trajectory)
PropagationResult calculateTrackParametersPropagateHelper(const EventContext &, const TrackParameters &, const GXFTrackState &, PropDirection, const MagneticFieldProperties &, bool, bool) const
Helper method that encapsulates calls to the propagator tool in the calculateTrackParameters() method...
ToolHandle< IUpdator > m_updator
static void fillBfromMeasurements(const Cache &cache, GXFTrajectory &trajectory, Amg::VectorX &b)
std::unique_ptr< const TrackParameters > makeTrackFindPerigeeParameters(const EventContext &, Cache &, GXFTrajectory &, const ParticleHypothesis) const
SG::ReadCondHandleKey< AtlasFieldCacheCondObj > m_field_cache_key
Gaudi::Property< bool > m_domeastrackpar
virtual int iterationsOfLastFit() const
Gaudi::Property< int > m_fixbrem
Track * backupCombinationStrategy(const EventContext &ctx, Cache &, const Track &, const Track &, GXFTrajectory &, std::vector< MaterialEffectsOnTrack > &) const
GlobalChi2Fitter(const std::string &, const std::string &, const IInterface *)
Track * mainCombinationStrategy(const EventContext &ctx, Cache &, const Track &, const Track &, GXFTrajectory &, std::vector< MaterialEffectsOnTrack > &) const
static void compensatePhiWeights(Cache &cache, GXFTrajectory &trajectory, Amg::SymMatrixX &a)
const TrackingGeometry * trackingGeometry(Cache &cache, const EventContext &ctx) const
Gaudi::Property< int > m_maxit
PropagationResult calculateTrackParametersPropagate(const EventContext &, const TrackParameters &, const GXFTrackState &, PropDirection, const MagneticFieldProperties &, bool, bool) const
Propagate onto a track state, collecting new track parameters, and optionally the Jacobian and possib...
ToolHandle< IExtrapolator > m_extrapolator
void makeProtoStateFromMeasurement(Cache &, GXFTrajectory &, const MeasurementBase *, const TrackParameters *trackpar=nullptr, bool isoutlier=false, int index=-1) const
std::optional< TransportJacobian > numericalDerivatives(const EventContext &ctx, const TrackParameters *, const Surface &, PropDirection, const MagneticFieldProperties &) const
void addIDMaterialFast(const EventContext &ctx, Cache &cache, GXFTrajectory &track, const TrackParameters *parameters, ParticleHypothesis part) const
A faster strategy for adding scatter material to tracks, works only for inner detector tracks.
static void fillAfromScatterers(GXFTrajectory &trajectory, Amg::SymMatrixX &a)
void updatePixelROTs(GXFTrajectory &, Amg::SymMatrixX &, Amg::VectorX &, const EventContext &evtctx) const
Update the Pixel ROT using the current trajectory/local track parameters.
Gaudi::Property< double > m_outlcut
Gaudi::Property< bool > m_calomat
Gaudi::Property< bool > m_fillderivmatrix
ToolHandle< IMaterialEffectsOnTrackProvider > m_calotoolparam
ToolHandle< IRIO_OnTrackCreator > m_ROTcreator
GXFTrajectory * runTrackCleanerSilicon(const EventContext &ctx, Cache &, GXFTrajectory &, Amg::SymMatrixX &, Amg::SymMatrixX &, Amg::VectorX &, bool) const
std::unique_ptr< const TrackParameters > makePerigee(Cache &, const TrackParameters &, const ParticleHypothesis) const
virtual std::unique_ptr< Track > fit(const EventContext &ctx, const PrepRawDataSet &, const TrackParameters &, const RunOutlierRemoval runOutlier=false, const ParticleHypothesis matEffects=nonInteracting) const override final
Gaudi::Property< double > m_chi2cut
void throwFailedToGetTrackingGeomtry() const
Gaudi::Property< bool > m_extensioncuts
virtual void setMinIterations(int)
bool ensureValidEntranceCalo(const EventContext &ctx, Cache &cache) const
static void addMaterialGetLayers(Cache &cache, std::vector< std::pair< const Layer *, const Layer * > > &layers, std::vector< std::pair< const Layer *, const Layer * > > &uplayers, const std::vector< std::unique_ptr< GXFTrackState > > &states, GXFTrackState &first, GXFTrackState &last, const TrackParameters *refpar, bool hasmat)
Collect all possible layers that a given track could have passed through.
static std::optional< std::pair< Amg::Vector3D, double > > addMaterialFindIntersectionCyl(Cache &cache, const CylinderSurface &surface, const TrackParameters ¶m1, const TrackParameters ¶m2, const ParticleHypothesis mat)
Find the intersection of a set of track parameters onto a cylindrical surface.
void runTrackCleanerTRT(Cache &, GXFTrajectory &, Amg::SymMatrixX &, Amg::VectorX &, Amg::SymMatrixX &, bool, bool, int, const EventContext &ctx) const
static void calculateDerivatives(GXFTrajectory &)
Gaudi::Property< bool > m_decomposesegments
void calculateTrackErrors(GXFTrajectory &, Amg::SymMatrixX &, bool) const
static bool correctAngles(double &, double &)
bool ensureValidEntranceMuonSpectrometer(const EventContext &ctx, Cache &cache) const
bool isMuonTrack(const Track &) const
std::vector< std::unique_ptr< TrackParameters > > holesearchExtrapolation(const EventContext &ctx, const TrackParameters &src, const GXFTrackState &dst, PropDirection propdir) const
Helper method which performs an extrapolation with additional logic for hole search.
virtual StatusCode finalize() override
Gaudi::Property< bool > m_useCaloTG
std::vector< std::reference_wrapper< GXFTrackState > > holeSearchStates(GXFTrajectory &trajectory) const
Extracts a collection of track states which are important for hole search.
Track * fitIm(const EventContext &ctx, Cache &cache, const Track &inputTrack, const RunOutlierRemoval runOutlier, const ParticleHypothesis matEffects) const
Gaudi::Property< bool > m_redoderivs
void initFieldCache(const EventContext &ctx, Cache &cache) const
Initialize a field cache inside a fit cache object.
Gaudi::Property< int > m_maxitPixelROT
ToolHandle< IMaterialEffectsOnTrackProvider > m_calotool
ToolHandle< IMaterialEffectsUpdator > m_matupdator
void holeSearchHelper(const std::vector< std::unique_ptr< TrackParameters > > &hc, std::set< Identifier > &id_set, std::set< Identifier > &sct_set, TrackHoleCount &rv, bool count_holes, bool count_dead) const
Helper method for the hole search that does the actual counting of holes and dead modules.
ToolHandle< IResidualPullCalculator > m_residualPullCalculator
void addMaterialUpdateTrajectory(Cache &cache, GXFTrajectory &track, int offset, std::vector< std::pair< const Layer *, const Layer * > > &layers, const TrackParameters *ref1, const TrackParameters *ref2, ParticleHypothesis mat) const
Given layer information, probe those layers for scatterers and add them to a track.
ToolHandle< IEnergyLossUpdator > m_elosstool
void fillDerivatives(GXFTrajectory &traj) const
Gaudi::Property< bool > m_numderiv
SG::ReadCondHandleKey< TrackingGeometry > m_trackingGeometryReadKey
Gaudi::Property< bool > m_acceleration
static bool tryToWeightAfromMaterial(Cache &cache, GXFTrajectory &trajectory, Amg::SymMatrixX &a, const bool doDeriv, const int it, const double oldRedChi2, const double newRedChi2)
const AtlasDetectorID * m_DetID
static void fillAfromMeasurements(const Cache &cache, GXFTrajectory &trajectory, Amg::SymMatrixX &a)
virtual StatusCode initialize() override
static void makeTrackFillDerivativeMatrix(Cache &, GXFTrajectory &)
Gaudi::Property< bool > m_holeSearch
FitterStatusCode runIteration(const EventContext &ctx, Cache &cache, GXFTrajectory &trajectory, const int it, Amg::SymMatrixX &a, Amg::VectorX &b, Amg::SymMatrixX &lu, bool &doDeriv) const
Gaudi::Property< bool > m_createSummary
ToolHandle< IMultipleScatteringUpdator > m_scattool
Interface class IPropagators It inherits from IAlgTool.
LayerIndex for the identification of layers in a simplified detector geometry of Cylinders and Discs.
int value() const
layerIndex expressed in an integer
double oppositePostFactor() const
Return method for post update material description of the Layer along normalvector.
double alongPreFactor() const
Return method for pre update material description of the Layer along normalvector.
double alongPostFactor() const
Return method for post update material description of the Layer along normalvector.
double oppositePreFactor() const
Return method for pre update material description of the Layer along normalvector.
Base Class for a Detector Layer in the Tracking realm.
virtual const Surface & surfaceRepresentation() const =0
Transforms the layer into a Surface representation for extrapolation.
const MaterialProperties * fullUpdateMaterialProperties(const TrackParameters &par) const
getting the MaterialProperties back - for full update
const LayerMaterialProperties * layerMaterialProperties() const
getting the LayerMaterialProperties including full/pre/post update
int parameterKey() const
Identifier key for matrix expansion/reduction.
bool contains(ParamDefs par) const
The simple check for the clients whether the parameter is contained.
magnetic field properties to steer the behavior of the extrapolation
base class to integrate material effects on Trk::Track in a flexible way.
const Surface & associatedSurface() const
returns the surface to which these m.eff. are associated.
double thicknessInX0() const
returns the actually traversed material .
@ MATERIAL_EFFECTS_ON_TRACK
virtual MaterialEffectsDerivedType derivedType() const =0
Returns the concrete derived type.
represents the full description of deflection and e-loss of a track in material.
const EnergyLoss * energyLoss() const
returns the energy loss object.
const ScatteringAngles * scatteringAngles() const
returns the MCS-angles object.
Material with information about thickness of material.
float thicknessInX0() const
Return the radiationlength fraction.
float thickness() const
Return the thickness in mm.
This class is the pure abstract base class for all fittable tracking measurements.
virtual MeasurementBase * clone() const =0
Pseudo-Constructor.
const LocalParameters & localParameters() const
Interface method to get the LocalParameters.
virtual const Surface & associatedSurface() const =0
Interface method to get the associated Surface.
virtual bool type(MeasurementBaseType::Type type) const =0
Interface method checking the type.
virtual const Amg::Vector3D & globalPosition() const =0
Interface method to get the global Position.
const Amg::MatrixX & localCovariance() const
Interface method to get the localError.
const Amg::Vector3D & momentum() const
Access method for the momentum.
virtual ParametersBase< DIM, T > * clone() const override=0
clone method for polymorphic deep copy
const Amg::Vector3D & position() const
Access method for the position.
virtual const Surface & associatedSurface() const override=0
Access to the Surface associated to the Parameters.
double pT() const
Access method for transverse momentum.
Class describing the Line to which the Perigee refers to.
virtual DistanceSolution straightLineDistanceEstimate(const Amg::Vector3D &pos, const Amg::Vector3D &dir) const override final
fast straight line distance evaluation to Surface
virtual Surface::ChargedTrackParametersUniquePtr createUniqueTrackParameters(double l1, double l2, double phi, double theta, double qop, std::optional< AmgSymMatrix(5)> cov=std::nullopt) const override final
Use the Surface as a ParametersBase constructor, from local parameters - charged.
Class for a planaer rectangular or trapezoidal surface in the ATLAS detector.
const Amg::Vector2D & localPosition() const
return the local position reference
virtual bool type(PrepRawDataType type) const
Interface method checking the type.
const Amg::MatrixX & localCovariance() const
return const ref to the error matrix
Class to handle RIO On Tracks ROT) for InDet and Muons, it inherits from the common MeasurementBase.
virtual const Surface & associatedSurface() const override=0
returns the surface for the local to global transformation
virtual const Trk::PrepRawData * prepRawData() const =0
returns the PrepRawData (also known as RIO) object to which this RIO_OnTrack is associated.
Identifier identify() const
return the identifier -extends MeasurementBase
virtual const Amg::Vector3D & globalPosition() const override=0
Interface method to get the global Position.
@ Biased
RP with track state including the hit.
represents a deflection of the track caused through multiple scattering in material.
Base class for all TrackSegment implementations, extends the common MeasurementBase.
const MeasurementBase * measurement(unsigned int) const
returns the Trk::MeasurementBase objects depending on the integer
unsigned int numberOfMeasurementBases() const
Return the number of contained Trk::MeasurementBase (s)
Class for a StraightLineSurface in the ATLAS detector to describe dirft tube and straw like detectors...
virtual double r() const =0
Interface method for the maximal extension or the radius.
virtual BoundsType type() const =0
Return the bounds type - for persistency optimization.
Abstract Base Class for tracking surfaces.
virtual DistanceSolution straightLineDistanceEstimate(const Amg::Vector3D &pos, const Amg::Vector3D &dir) const =0
fast straight line distance evaluation to Surface
virtual ChargedTrackParametersUniquePtr createUniqueTrackParameters(double l1, double l2, double phi, double theat, double qop, std::optional< AmgSymMatrix(5)> cov=std::nullopt) const =0
Use the Surface as a ParametersBase constructor, from local parameters - charged.
virtual const Amg::Vector3D & normal() const
Returns the normal vector of the Surface (i.e.
const Amg::Transform3D & transform() const
Returns HepGeom::Transform3D by reference.
Identifier associatedDetectorElementIdentifier() const
return Identifier of the associated Detector Element
virtual constexpr SurfaceType type() const =0
Returns the Surface type to avoid dynamic casts.
virtual const SurfaceBounds & bounds() const =0
Surface Bounds method.
const Trk::Layer * associatedLayer() const
return the associated Layer
const Amg::Vector3D & center() const
Returns the center position of the Surface.
Contains information about the 'fitter' of this track.
bool trackProperties(const TrackProperties &property) const
Access methods for track properties.
@ GlobalChi2Fitter
Track's from Thijs' global chi^2 fitter.
@ Unknown
Track fitter not defined.
ParticleHypothesis particleHypothesis() const
Returns the particle hypothesis used for Track fitting.
const TrackFitter & trackFitter() const
Access methods for track fitter.
@ BremFit
A brem fit was performed on this track.
@ StraightTrack
A straight track.
@ BremFitSuccessful
A brem fit was performed on this track and this fit was successful.
@ SlimmedTrack
A slimmed track.
represents the track state (measurement, material, fit parameters and quality) at a surface.
const MeasurementBase * measurementOnTrack() const
returns MeasurementBase const overload
const TrackParameters * trackParameters() const
return ptr to trackparameters const overload
bool type(const TrackStateOnSurfaceType type) const
Use this method to find out if the TSoS is of a certain type: i.e.
@ Measurement
This is a measurement, and will at least contain a Trk::MeasurementBase.
@ Perigee
This represents a perigee, and so will contain a Perigee object only.
@ BremPoint
This represents a brem point on the track, and so will contain TrackParameters and MaterialEffectsBas...
@ Outlier
This TSoS contains an outlier, that is, it contains a MeasurementBase/RIO_OnTrack which was not used ...
@ InertMaterial
This represents inert material, and so will contain MaterialEffectsBase.
@ Scatterer
This represents a scattering point on the track, and so will contain TrackParameters and MaterialEffe...
@ CaloDeposit
This TSOS contains a CaloEnergy object.
const MaterialEffectsBase * materialEffectsOnTrack() const
return material effects const overload
const Trk::Surface & surface() const
return associated surface
const Trk::TrackStates * trackStateOnSurfaces() const
return a pointer to a const DataVector of const TrackStateOnSurfaces.
const DataVector< const MeasurementBase > * measurementsOnTrack() const
return a pointer to a vector of MeasurementBase (NOT including any that come from outliers).
const DataVector< const TrackParameters > * trackParameters() const
Return a pointer to a vector of TrackParameters.
const TrackInfo & info() const
Returns a const ref to info of a const tracks.
const Perigee * perigeeParameters() const
return Perigee.
const FitQuality * fitQuality() const
return a pointer to the fit quality const-overload
The TrackingGeometry class is the owner of the constructed TrackingVolumes.
Full Volume description used in Tracking, it inherits from Volume to get the geometrical structure,...
const LayerArray * confinedLayers() const
Return the subLayer array.
LayerIntersection< Amg::Vector3D > closestMaterialLayer(const Amg::Vector3D &gp, const Amg::Vector3D &dir, PropDirection pDir=alongMomentum, const BoundaryCheck &bchk=true) const
Return the closest layer with material description.
const TrackingVolumeArray * confinedVolumes() const
Return the subLayer array.
std::vector< std::shared_ptr< BoundarySurface< TrackingVolume > > > & boundarySurfaces()
Method to return the BoundarySurfaces.
This is the base class for all tracking detector elements with read-out relevant information.
virtual Identifier identify() const =0
Identifier.
bool inside(const Amg::Vector3D &gp, double tol=0.) const
Inside() method for checks.
double chi2(TH1 *h0, TH1 *h1)
std::string head(std::string s, const std::string &pattern)
head of a string
Eigen::Matrix< double, 3, 3 > RotationMatrix3D
Eigen::Matrix< double, Eigen::Dynamic, Eigen::Dynamic > MatrixX
Dynamic Matrix - dynamic allocation.
Eigen::Matrix< double, Eigen::Dynamic, Eigen::Dynamic > SymMatrixX
Eigen::Affine3d Transform3D
Eigen::Matrix< double, 3, 1 > Vector3D
Eigen::Matrix< double, Eigen::Dynamic, 1 > VectorX
Dynamic Vector - dynamic allocation.
double getDistance(const xAOD::Vertex *vtx1, const xAOD::Vertex *vtx2)
@ PseudoMeasurementOnTrack
constexpr double mass[PARTICLEHYPOTHESES]
the array of masses
MeasurementType
enum describing the flavour of MeasurementBase
Ensure that the ATLAS eigen extensions are properly loaded.
PropDirection
PropDirection, enum for direction of the propagation.
std::vector< const MeasurementBase * > MeasurementSet
vector of fittable measurements
ParametersT< TrackParametersDim, Charged, PerigeeSurface > Perigee
ParametersT< TrackParametersDim, Charged, StraightLineSurface > AtaStraightLine
@ DeadElement
outside the element
bool consistentSurfaces(U)
std::unique_ptr< T > unique_clone(const T *v)
bool RunOutlierRemoval
switch to toggle quality processing after fit
@ NoField
Field is set to 0., 0., 0.,.
@ FullField
Field is set to be realistic, but within a given Volume.
@ z
global position (cartesian)
@ u
Enums for curvilinear frames.
@ loc2
generic first and second local coordinate
std::pair< double, ParamDefs > DefinedParameter
Typedef to of a std::pair<double, ParamDefs> to identify a passed-through double as a specific type o...
ParticleHypothesis
Enumeration for Particle hypothesis respecting the interaction with material.
std::pair< long int, long int > indices
BinnedArray< TrackingVolume > TrackingVolumeArray
simply for the eye
ParametersBase< TrackParametersDim, Charged > TrackParameters
ParametersT< TrackParametersDim, Charged, PlaneSurface > AtaPlane
std::vector< const PrepRawData * > PrepRawDataSet
vector of clusters and drift circles
AmgSymMatrix(5) &GXFTrackState
@ numberOfSCTHoles
number of Holes in both sides of a SCT module
@ numberOfSCTDeadSensors
number of TRT hits
@ numberOfPixelHoles
number of pixels which have a ganged ambiguity.
@ numberOfPixelDeadSensors
number of pixel hits with broad errors (width/sqrt(12))
void sort(typename DataModel_detail::iterator< DVL > beg, typename DataModel_detail::iterator< DVL > end)
Specialization of sort for DataVector/List.
void stable_sort(DataModel_detail::iterator< DVL > beg, DataModel_detail::iterator< DVL > end)
Specialization of stable_sort for DataVector/List.
void reverse(typename DataModel_detail::iterator< DVL > beg, typename DataModel_detail::iterator< DVL > end)
Specialization of reverse for DataVector/List.
double deltaPhi(double phiA, double phiB)
delta Phi in range [-pi,pi[
std::vector< std::unique_ptr< const std::vector< const TrackStateOnSurface * >, void(*)(const std::vector< const TrackStateOnSurface * > *) > > m_matTempStore
void incrementFitStatus(enum FitterStatusType status)
Amg::SymMatrixX m_fullcovmat
bool m_getmaterialfromtrack
std::vector< const Trk::Layer * > m_posdiscs
std::vector< int > m_lastmeasurement
std::vector< const Trk::Layer * > m_negdiscs
std::vector< const Trk::Layer * > m_barrelcylinders
MagField::AtlasFieldCache m_field_cache
static void objVectorDeleter(const std::vector< const T * > *ptr)
std::vector< int > m_firstmeasurement
std::array< unsigned int, S_MAX_VALUE > m_fit_status
std::vector< double > m_phiweight
const TrackingVolume * m_msEntrance
const TrackingVolume * m_caloEntrance
FitterStatusCode m_fittercode
static constexpr std::array< ParamDefs, 6 > pardef
Constructor.