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 delete alignCache.m_derivMatrix;
1830 alignCache.m_derivMatrix =
nullptr;
1832 delete alignCache.m_fullCovarianceMatrix;
1833 alignCache.m_fullCovarianceMatrix =
nullptr;
1834 alignCache.m_iterationsOfLastFit = 0;
1837 fitIm(ctx, cache, inputTrack, runOutlier, matEffects);
1838 if(newTrack !=
nullptr){
1843 alignCache.m_iterationsOfLastFit = cache.
m_lastiter;
1850 const EventContext& ctx,
1852 const Track& inputTrack,
1857 ATH_MSG_DEBUG(
"--> entering GlobalChi2Fitter::fit(Track,,)");
1873 ATH_MSG_WARNING(
"Track without track parameters, cannot perform fit");
1881 if (minpar ==
nullptr) {
1899 const Surface *firsthitsurf =
nullptr;
1900 const Surface *lasthitsurf =
nullptr;
1902 bool hasmuon =
false;
1903 bool iscombined =
false;
1904 bool seenphimeas =
false;
1908 for (; itStates != endState; ++itStates) {
1909 const auto *
const pMeasurement = (**itStates).measurementOnTrack();
1911 (pMeasurement ==
nullptr) &&
1912 ((**itStates).materialEffectsOnTrack() !=
nullptr) &&
1918 if (pMeasurement !=
nullptr) {
1919 const Surface *surf = &pMeasurement->associatedSurface();
1928 if (firsthitsurf ==
nullptr) {
1929 firsthitsurf = surf;
1932 if (
m_DetID->is_indet(hitid)) {
1937 if ((**itStates).trackParameters() !=
nullptr) {
1938 lastidpar = (**itStates).trackParameters();
1939 if (firstidpar ==
nullptr) {
1940 firstidpar = lastidpar;
1946 const double dotprod1 = measdir.dot(
Amg::Vector3D(0, 0, 1));
1948 if (std::abs(dotprod1) < .5 && std::abs(dotprod2) < .5) {
1950 if (std::abs(surf->
center().z()) > 13000) {
1953 if (surf->
center().perp() > 9000 && std::abs(surf->
center().z()) < 13000) {
1965 if (iscombined && seenphimeas && (phiem || phibo)) {
1981 if (firstidpar == lastidpar) {
1982 firstidpar = lastidpar =
nullptr;
1992 (firstidpar !=
nullptr)
2003 const bool tmpsirecal = cache.
m_sirecal;
2004 std::unique_ptr<Track> tmptrack =
nullptr;
2015 tmptrack.reset(
myfit(ctx, cache, trajectory, *minpar,
false, matEffects));
2018 if (tmptrack ==
nullptr) {
2027 bool isbrem =
false;
2029 unsigned int n_brem=0;
2031 for (std::unique_ptr<GXFTrackState> & state : trajectory.
trackStates()) {
2034 if (meff !=
nullptr) {
2040 const double p = 1. / std::abs(layerpars->parameters()[
Trk::qOverP] - .0005 * meff->
delta_p());
2042 std::optional<Amg::Vector2D> locpos(state->associatedSurface().globalToLocal(layerpars->
position()));
2043 const Amg::Vector3D layerNormal(state->associatedSurface().normal(*locpos));
2044 double costracksurf = 1.;
2046 costracksurf = std::abs(layerNormal.dot(layerpars->
momentum().unit()));
2048 const double oldde = meff->
deltaE();
2050 std::unique_ptr<EnergyLoss> eloss;
2051 double sigmascat = 0;
2053 if (matprop !=
nullptr) {
2054 eloss = std::make_unique<EnergyLoss>(
2055 m_elosstool->energyLoss(*matprop, p, 1. / costracksurf,
2057 sigmascat = std::sqrt(
m_scattool->sigmaSquare(*matprop, p, 1. / costracksurf, matEffects));
2059 if (eloss !=
nullptr) {
2064 sigmascat = std::sqrt(
m_scattool->sigmaSquare(tmpprop, p, 1. / costracksurf, matEffects));
2068 sigmascat / std::sin(layerpars->parameters()[
Trk::theta]),
2082 }
else if (eloss !=
nullptr) {
2093 refpars[0], refpars[1], refpars[2], refpars[3], refpars[4], std::nullopt
2101 trajectory.
brems().clear();
2103 trajectory.
brems().resize(1);
2104 trajectory.
brems()[0] = bremdp;
2115 std::unique_ptr<Track> track(
myfit(ctx, cache, trajectory, *minpar, runOutlier, matEffects));
2117 bool pseudoupdated =
false;
2119 if ((track !=
nullptr) && hasid && hasmuon) {
2120 for (std::unique_ptr<GXFTrackState> & pseudostate : trajectory.
trackStates()) {
2122 (pseudostate ==
nullptr) ||
2124 pseudostate->fitQuality().chiSquared() < 10
2130 const std::unique_ptr<const TrackParameters> updpar(
m_updator->removeFromState(
2132 pseudostate->measurement()->localParameters(),
2133 pseudostate->measurement()->localCovariance()
2136 if (updpar ==
nullptr) {
2141 covMatrix(0, 0) = 100;
2143 std::unique_ptr<const PseudoMeasurementOnTrack> newpseudo = std::make_unique<const PseudoMeasurementOnTrack>(
2145 std::move(covMatrix),
2149 pseudostate->setMeasurement(std::move(newpseudo));
2151 errors[0] = errors[2] = errors[3] = errors[4] = -1;
2153 pseudostate->setMeasurementErrors(errors);
2154 pseudoupdated =
true;
2157 if (pseudoupdated) {
2160 track.reset(
myfit(ctx, cache, trajectory, *track->perigeeParameters(),
false,
muon));
2170 if (track !=
nullptr) {
2173 track->info().addPatternReco(old_info);
2176 return track.release();
2179 std::unique_ptr<Track>
2181 const EventContext& ctx,
2187 ATH_MSG_DEBUG(
"--> entering GlobalChi2Fitter::fit(PRDS,TP,)");
2190 for (
const auto *prd : prds) {
2191 const Surface & prdsurf = (*prd).detectorElement()->surface((*prd).identify());
2196 plsurf =
static_cast < const PlaneSurface *
>(&prdsurf);
2203 if ((slsurf ==
nullptr) && (plsurf ==
nullptr)) {
2204 ATH_MSG_ERROR(
"Surface is neither PlaneSurface nor StraightLineSurface!");
2209 }
else if (slsurf !=
nullptr) {
2218 }
else if (plsurf !=
nullptr) {
2219 if (param.covariance() !=
nullptr) {
2241 if (rot !=
nullptr) {
2242 rots.push_back(rot);
2246 std::unique_ptr<Track> track =
2247 fit(ctx, rots, param, runOutlier, matEffects);
2249 for (
const auto *rot : rots) {
2256 std::unique_ptr<Track>
2258 const EventContext& ctx,
2259 const Track& inputTrack,
2264 ATH_MSG_DEBUG(
"--> entering GlobalChi2Fitter::fit(Track,Meas'BaseSet,,)");
2279 if (minpar ==
nullptr) {
2297 for (; itStates != endState; ++itStates) {
2302 (trajectory.
trackStates().back()->materialEffects() !=
nullptr) &&
2303 trajectory.
trackStates().back()->materialEffects()->sigmaDeltaE() > 50.001
2305 trajectory.
trackStates().back()->materialEffects()->setKink(
true);
2311 for (
const auto & measBase : addMeasColl) {
2312 if (measBase ==
nullptr) {
2313 ATH_MSG_WARNING(
"There is an empty MeasurementBase object in the track! Skip this object..");
2321 std::unique_ptr<Track> track(
myfit(ctx, cache, trajectory, *minpar, runOutlier, matEffects));
2324 if (track !=
nullptr) {
2325 const double oldqual =
2330 const double newqual =
2331 track->fitQuality()->numberDoF() != 0 ?
2332 track->fitQuality()->chiSquared() / track->fitQuality()->numberDoF() :
2335 if (
m_extensioncuts && runOutlier && newqual > 2 && newqual > 2 * oldqual) {
2339 track.reset(
nullptr);
2343 if (track !=
nullptr) {
2352 std::unique_ptr<Track>
2354 const EventContext& ctx,
2360 ATH_MSG_DEBUG(
"--> entering GlobalChi2Fitter::fit(Track,PRDS,)");
2364 for (
const auto *prd : prds) {
2365 const Surface & prdsurf = (*prd).detectorElement()->surface((*prd).identify());
2368 std::unique_ptr<const TrackParameters>
const trackparForCorrect(
2384 rot =
m_ROTcreator->correct(*prd, *trackparForCorrect, ctx);
2387 if (rot !=
nullptr) {
2388 rots.push_back(rot);
2392 std::unique_ptr<Track> track =
fit(ctx,intrk, rots, runOutlier, matEffects);
2394 for (
const auto *rot : rots) {
2402 const EventContext& ctx,
2408 ATH_MSG_DEBUG(
"--> entering GlobalChi2Fitter::fit(Meas'BaseSet,,)");
2421 for (
const auto *itSet : rots) {
2422 if (itSet ==
nullptr) {
2423 ATH_MSG_WARNING(
"There is an empty MeasurementBase object in the track! Skip this object..");
2429 std::unique_ptr<const TrackParameters> startpar(param.
clone());
2432 matEffects ==
muon &&
2438 myfit(ctx,cache, trajectory, *startpar, runOutlier, matEffects);
2454 covMatrix(0, 0) = 100;
2456 std::unique_ptr<const PseudoMeasurementOnTrack> newpseudo = std::make_unique<const PseudoMeasurementOnTrack>(
2458 std::move(covMatrix),
2462 trajectory.
trackStates().front()->setMeasurement(std::move(newpseudo));
2467 covMatrix(0, 0) = 100;
2469 std::unique_ptr<const PseudoMeasurementOnTrack> newpseudo = std::make_unique<const PseudoMeasurementOnTrack>(
2471 std::move(covMatrix),
2475 trajectory.
trackStates().back()->setMeasurement(std::move(newpseudo));
2482 refpars[0], refpars[1], refpars[2], refpars[3], refpars[4], std::nullopt
2487 myfit(ctx,cache, trajectory, *startpar, runOutlier, matEffects);
2498 refpars[0], refpars[1], refpars[2], refpars[3], refpars[4], std::nullopt
2505 firstpar = trajectory.
trackStates().front()->trackParameters();
2508 covMatrix(0, 0) = 100;
2510 std::unique_ptr<const PseudoMeasurementOnTrack> newpseudo = std::make_unique<const PseudoMeasurementOnTrack>(
2512 std::move(covMatrix),
2516 trajectory.
trackStates().front()->setMeasurement(std::move(newpseudo));
2518 errors[0] = errors[2] = errors[3] = errors[4] = -1;
2520 trajectory.
trackStates().front()->setMeasurementErrors(errors);
2524 lastpar = trajectory.
trackStates().back()->trackParameters();
2527 covMatrix(0, 0) = 100;
2529 std::unique_ptr<const PseudoMeasurementOnTrack> newpseudo = std::make_unique<const PseudoMeasurementOnTrack>(
2531 std::move(covMatrix),
2535 trajectory.
trackStates().back()->setMeasurement(std::move(newpseudo));
2537 errors[0] = errors[2] = errors[3] = errors[4] = -1;
2539 trajectory.
trackStates().back()->setMeasurementErrors(errors);
2543 std::unique_ptr<Track> track;
2545 if (startpar !=
nullptr) {
2546 track.reset(
myfit(ctx,cache, trajectory, *startpar, runOutlier, matEffects));
2549 if (track !=
nullptr) {
2579 std::unique_ptr<GXFMaterialEffects> newmeff;
2587 newmeff = std::make_unique<GXFMaterialEffects>(*meff);
2591 const double sigmascat = std::sqrt(
m_scattool->sigmaSquare(
2612 newmeff = std::make_unique<GXFMaterialEffects>(newmeot);
2628 ((**trajectory.
trackStates().rbegin()).trackParameters() !=
nullptr)
2630 const double delta_p = 1000 * (
2632 (**trajectory.
trackStates().rbegin()).trackParameters()->
2636 newmeff->setdelta_p(delta_p);
2672 seg =
static_cast<const Segment *
>(measbase);
2682 for (
int i = 0; i <
imax; i++) {
2685 std::unique_ptr<GXFTrackState> ptsos = std::make_unique<GXFTrackState>(
2686 std::unique_ptr<const MeasurementBase>(measbase2->
clone()),
2687 std::unique_ptr<const TrackParameters>(newtrackpar !=
nullptr ? newtrackpar->
clone() :
nullptr)
2690 double sinstereo = 0;
2692 errors[0] = errors[1] = errors[2] = errors[3] = errors[4] = -1;
2703 bool measphi =
false;
2706 bool rotated =
false;
2709 if (
m_DetID->is_pixel(hitid)) {
2711 }
else if (
m_DetID->is_sct(hitid)) {
2712 if (covmat.cols() != 1 && covmat(1, 0) != 0) {
2716 }
else if (
m_DetID->is_trt(hitid)) {
2726 }
else if (
m_DetID->is_mdt(hitid)) {
2728 }
else if (
m_DetID->is_tgc(hitid)) {
2733 }
else if (
m_DetID->is_csc(hitid)) {
2735 }
else if (
m_DetID->is_mm(hitid)) {
2737 }
else if (
m_DetID->is_stgc(hitid)) {
2743 const auto [covEigenValueSmall, covStereoAngle] = principalComponentAnalysis2x2(covmat);
2744 errors[0] = std::sqrt(covEigenValueSmall);
2745 sinstereo = std::sin(covStereoAngle);
2747 errors[0] = std::sqrt(covmat(0, 0));
2749 errors[1] = std::sqrt(covmat(1, 1));
2760 const double dotprod1 = measdir.dot(
Amg::Vector3D(0, 0, 1));
2762 if (std::abs(dotprod1) < .5 && std::abs(dotprod2) < .5) {
2770 int param_index = 0;
2772 errors[0] = std::sqrt(covmat(0, 0));
2777 errors[1] = std::sqrt(covmat(param_index, param_index));
2782 errors[2] = std::sqrt(covmat(param_index, param_index));
2787 errors[3] = std::sqrt(covmat(param_index, param_index));
2792 errors[4] = std::sqrt(covmat(param_index, param_index));
2813 ptsos->setMeasurementErrors(errors);
2814 ptsos->setSinStereo(sinstereo);
2815 ptsos->setMeasurementType(hittype);
2816 ptsos->setMeasuresPhi(measphi);
2837 if (tvol ==
nullptr) {
2841 const Trk::BinnedArray < Trk::Layer > *confinedLayers = tvol->
confinedLayers();
2844 if (confinedLayers !=
nullptr) {
2846 for (
const auto & layer : confinedLayers->
arrayObjects()) {
2848 if (layer !=
nullptr) {
2853 if ((layIndex.
value() == 0) || (layer->layerMaterialProperties() ==
nullptr)) {
2864 disclay =
static_cast<const DiscLayer *
>(layer);
2867 if (disclay !=
nullptr) {
2868 if (disclay->
center().z() < 0) {
2873 }
else if (cyllay !=
nullptr) {
2884 for (
size_t ib = 0 ; ib < bsurf.size(); ++ib) {
2885 const Layer *layer = bsurf[ib]->surfaceRepresentation().materialLayer();
2887 if (layer ==
nullptr)
continue;
2891 if ((layIndex.
value() == 0) || (layer->layerMaterialProperties() ==
nullptr)) {
2898 cylsurf =
static_cast<const CylinderSurface *
>(&layer->surfaceRepresentation());
2903 discsurf =
static_cast<const DiscSurface *
>(&layer->surfaceRepresentation());
2905 if (discsurf !=
nullptr) {
2907 discsurf->
center().z() < 0 &&
2912 discsurf->
center().z() > 0 &&
2918 (cylsurf !=
nullptr) &&
2924 if ((cylsurf ==
nullptr) && (discsurf ==
nullptr)) {
2931 if (confinedVolumes !=
nullptr) {
2932 for (
const auto & volume : confinedVolumes->
arrayObjects()) {
2933 if (volume !=
nullptr) {
2949 std::optional<std::pair<Amg::Vector3D, double>>
2962 const double * pos = parforextrap.
position().data();
2965 const double sinphi = std::sin(parforextrap.parameters()[
Trk::phi0]);
2966 const double cosphi = std::cos(parforextrap.parameters()[
Trk::phi0]);
2967 const double sintheta = std::sin(parforextrap.parameters()[
Trk::theta]);
2968 const double costheta = std::cos(parforextrap.parameters()[
Trk::theta]);
2972 const double r = (std::abs(currentqoverp) > 1e-10) ? -sintheta / (currentqoverp * 300. * field[2]) : 1e6;
2973 const double xc = parforextrap.
position().x() -
r * sinphi;
2974 const double yc = parforextrap.
position().y() +
r * cosphi;
2975 const double phi0 = std::atan2(parforextrap.
position().y() - yc, parforextrap.
position().x() - xc);
2976 const double z0 = parforextrap.
position().z();
2977 const double delta_s = (surf.
center().
z() -
z0) / costheta;
2978 const double delta_phi = delta_s * sintheta /
r;
2979 const double x = xc + std::abs(
r) * std::cos(
phi0 + delta_phi);
2980 const double y = yc + std::abs(
r) * std::sin(
phi0 + delta_phi);
2982 const double perp = intersect.perp();
2989 const double costracksurf = std::abs(costheta);
2991 return std::make_pair(intersect, costracksurf);
2994 std::optional<std::pair<Amg::Vector3D, double>>
3009 const double * pos = parforextrap.
position().data();
3012 const double sinphi = std::sin(parforextrap.parameters()[
Trk::phi0]);
3013 const double cosphi = std::cos(parforextrap.parameters()[
Trk::phi0]);
3014 const double sintheta = std::sin(parforextrap.parameters()[
Trk::theta]);
3015 const double costheta = std::cos(parforextrap.parameters()[
Trk::theta]);
3016 const double tantheta = std::tan(parforextrap.parameters()[
Trk::theta]);
3017 const double r = (std::abs(currentqoverp) > 1e-10) ? -sintheta / (currentqoverp * 300. * field[2]) : 1e6;
3018 const double xc = parforextrap.
position().x() -
r * sinphi;
3019 const double yc = parforextrap.
position().y() +
r * cosphi;
3020 const double phi0 = std::atan2(parforextrap.
position().y() - yc, parforextrap.
position().x() - xc);
3021 const double z0 = parforextrap.
position().z();
3022 const double d = xc * xc + yc * yc;
3023 const double rcyl = surf.
bounds().
r();
3024 double mysqrt = ((
r + rcyl) * (
r + rcyl) - d) * (d - (
r - rcyl) * (
r - rcyl));
3030 mysqrt = std::sqrt(mysqrt);
3031 double firstterm = xc / 2 + (xc * (rcyl * rcyl -
r *
r)) / (2 * d);
3032 double secondterm = (mysqrt * yc) / (2 * d);
3033 const double x1 = firstterm + secondterm;
3034 const double x2 = firstterm - secondterm;
3035 firstterm = yc / 2 + (yc * (rcyl * rcyl -
r *
r)) / (2 * d);
3036 secondterm = (mysqrt * xc) / (2 * d);
3037 const double y1 = firstterm - secondterm;
3038 const double y2 = firstterm + secondterm;
3041 const double dist1 = (
x - x1) * (
x - x1) + (
y - y1) * (
y - y1);
3042 const double dist2 = (
x - x2) * (
x - x2) + (
y - y2) * (
y - y2);
3044 if (dist1 < dist2) {
3052 const double phi1 = std::atan2(
y - yc,
x - xc);
3055 const double delta_z =
r * deltaphi / tantheta;
3056 const double z =
z0 + delta_z;
3067 const Amg::Vector3D trackdir(cos(phidir) * sintheta, std::sin(phidir) * sintheta, costheta);
3069 const double costracksurf = std::abs(normal.unit().dot(trackdir));
3071 return std::make_pair(intersect, costracksurf);
3078 std::vector<std::pair<const Layer *, const Layer *>> &layers,
3083 std::vector<std::unique_ptr<GXFTrackState>> & states = trajectory.
trackStates();
3084 std::vector<std::unique_ptr<GXFTrackState>> oldstates = std::move(states);
3087 states.reserve(oldstates.size() + layers.size());
3095 for (
int i = 0; i <= indexoffset; i++) {
3105 for (
int i = indexoffset + 1; i < (int) oldstates.size(); i++) {
3106 const double rmeas = oldstates[i]->position().perp();
3107 const double zmeas = oldstates[i]->position().z();
3115 while (layerindex < (
int) layers.size()) {
3117 double costracksurf = 0.0;
3118 const Layer *layer =
nullptr;
3127 if (layers[layerindex].first !=
nullptr) {
3131 layer = layers[layerindex].first;
3138 if (oldstates[i]->trackParameters() !=
nullptr) {
3139 const double rlayer = cylsurf->
bounds().
r();
3140 if (std::abs(rmeas - rlayer) < std::abs(parforextrap->
position().perp() - rlayer)) {
3141 parforextrap = oldstates[i]->trackParameters();
3151 std::tie(intersect, costracksurf) =
res.value();
3157 if (cylsurf->
bounds().
r() > rmeas)
break;
3158 }
else if (layers[layerindex].second !=
nullptr) {
3164 layer = layers[layerindex].second;
3167 if (oldstates[i]->trackParameters() !=
nullptr) {
3168 const double zlayer = discsurf->
center().z();
3169 if (std::abs(zmeas - zlayer) < std::abs(parforextrap->
position().z() - zlayer)) {
3170 parforextrap = oldstates[i]->trackParameters();
3175 std::tie(intersect, costracksurf) =
res.value();
3181 if (std::abs(discsurf->
center().z()) > std::abs(zmeas))
break;
3183 throw std::logic_error(
"Unhandled surface.");
3190 const MaterialProperties *matprop = layer->layerMaterialProperties()->fullMaterial(intersect);
3191 if (matprop ==
nullptr) {
3204 const double actualx0 = X0 / costracksurf;
3205 const double de = -std::abs(
3209 (
m_p != 0.0 ? std::abs(
m_p) : std::abs(1. / currentqoverp)),
3211 const double sintheta = std::sin(parforextrap->parameters()[
Trk::theta]);
3212 const double sigmascat = std::sqrt(
m_scattool->sigmaSquare(
3214 (
m_p != 0.0 ? std::abs(
m_p) : std::abs(1. / currentqoverp)),
3218 std::unique_ptr<GXFMaterialEffects> meff = std::make_unique<GXFMaterialEffects>();
3219 meff->setDeltaE(de);
3220 meff->setScatteringSigmas(sigmascat / sintheta, sigmascat);
3221 meff->setX0(actualx0);
3222 meff->setSurface(&layer->surfaceRepresentation());
3223 meff->setMaterialProperties(matprop);
3229 std::unique_ptr<EnergyLoss> eloss;
3232 eloss = std::make_unique<EnergyLoss>(
m_elosstool->energyLoss(
3234 (
m_p != 0.0 ? std::abs(
m_p) : std::abs(1. / currentqoverp)),
3239 if (eloss !=
nullptr) {
3240 meff->setSigmaDeltaE(eloss->sigmaDeltaE());
3245 meff->setDeltaE(-5);
3247 meff->setScatteringSigmas(0, 0);
3250 meff->setSigmaDeltaE(50);
3251 if (eloss !=
nullptr) {
3252 meff->setSigmaDeltaEPos(eloss->sigmaPlusDeltaE());
3253 meff->setSigmaDeltaENeg(eloss->sigmaMinusDeltaE());
3258 "X0: " << meff->x0() <<
" qoverp: " << currentqoverp <<
3259 " sigmascat " << meff->sigmaDeltaTheta() <<
" eloss: " << meff->deltaE() <<
3260 " sigma eloss: " << meff->sigmaDeltaE()
3267 std::unique_ptr<GXFTrackState> matstate = std::make_unique<GXFTrackState>(
3269 std::unique_ptr<const TrackParameters>()
3271 matstate->setPosition(intersect);
3287 std::vector<std::pair<const Layer *, const Layer *>> & layers,
3288 std::vector<std::pair<const Layer *, const Layer *>> & upstreamlayers,
3289 const std::vector<std::unique_ptr<GXFTrackState>> & oldstates,
3298 upstreamlayers.reserve(5);
3317 const double lastz = laststate->
position().z();
3318 const double lastr = laststate->
position().perp();
3324 const double tantheta = std::tan(refpar->parameters()[
Trk::theta]);
3325 const double slope = (tantheta != 0) ? 1 / tantheta : 0;
3331 std::vector < const Layer *>::const_iterator it;
3332 std::vector < const Layer *>::const_iterator itend;
3350 for (; it != itend; ++it) {
3355 if (std::abs((*it)->surfaceRepresentation().center().z()) > std::abs(lastz)) {
3363 const DiscBounds *discbounds =
static_cast<const DiscBounds *
> (&(*it)->surfaceRepresentation().bounds());
3368 if (discbounds->
rMax() < firstr || discbounds->
rMin() > lastr) {
3372 const double rintersect = firstr + ((*it)->surfaceRepresentation().center().
z() - firstz) / slope;
3375 rintersect < discbounds->rMin() - 50 ||
3376 rintersect > discbounds->
rMax() + 50
3386 if ((*it) == endlayer) {
3398 std::abs((*it)->surfaceRepresentation().center().z()) < std::abs(firstz) ||
3401 upstreamlayers.emplace_back((
Layer *)
nullptr, (*it));
3408 (*it) != startlayer &&
3409 (std::abs((*it)->surfaceRepresentation().center().z()) > std::abs(firstz2) ||
3410 (*it) == startlayer2)
3412 layers.emplace_back((
Layer *)
nullptr, (*it));
3424 if (barrelcylinder->surfaceRepresentation().bounds().r() > lastr) {
3431 const double zintersect = firstz + (barrelcylinder->surfaceRepresentation().bounds().
r() - firstr) * slope;
3433 if (std::abs(zintersect - barrelcylinder->surfaceRepresentation().center().z()) >
3434 ((
const CylinderSurface*)(&barrelcylinder->surfaceRepresentation()))->bounds().halflengthZ() + 50) {
3438 if (barrelcylinder == endlayer) {
3445 if (barrelcylinder->surfaceRepresentation().bounds().r() < firstr ||
3446 barrelcylinder == startlayer) {
3447 upstreamlayers.emplace_back(barrelcylinder, (
Layer*)
nullptr);
3450 if (barrelcylinder != startlayer &&
3451 (barrelcylinder->surfaceRepresentation().bounds().r() > firstr2 ||
3452 barrelcylinder == startlayer2)) {
3453 layers.emplace_back(barrelcylinder, (
Layer*)
nullptr);
3466 const EventContext& ctx,
3477 if (!caloEntranceIsValid) {
3504 addMaterial(ctx, cache, trajectory, refpar2, matEffects);
3520 bool hasmat =
false;
3521 int indexoffset = 0;
3522 int lastmatindex = 0;
3523 std::vector<std::unique_ptr<GXFTrackState>> & oldstates = trajectory.
trackStates();
3537 for (
int i = 0; i < (int) oldstates.size(); i++) {
3538 if (oldstates[i]->materialEffects() !=
nullptr) {
3547 if (firstsistate ==
nullptr) {
3548 if (oldstates[i]->trackParameters() ==
nullptr) {
3549 std::unique_ptr<const TrackParameters> tmppar(
m_propagator->propagateParameters(
3552 oldstates[i]->associatedSurface(),
3559 if (tmppar ==
nullptr)
return;
3561 oldstates[i]->setTrackParameters(std::move(tmppar));
3563 firstsistate = oldstates[i].get();
3565 lastsistate = oldstates[i].get();
3573 if (lastsistate ==
nullptr) {
3574 throw std::logic_error(
"No track state");
3584 std::unique_ptr<const TrackParameters> tmppar(
m_propagator->propagateParameters(
3593 if (tmppar ==
nullptr)
return;
3605 indexoffset = lastmatindex;
3620 std::vector<std::pair<const Layer *, const Layer *>> layers;
3621 std::vector<std::pair<const Layer *, const Layer *>> & upstreamlayers = trajectory.
upstreamMaterialLayers();
3626 addMaterialGetLayers(cache, layers, upstreamlayers, oldstates, *firstsistate, *lastsistate, refpar, hasmat);
3635 const EventContext& ctx,
3641 if (refpar2 ==
nullptr) {
3651 std::vector<std::unique_ptr<GXFTrackState>> & states = trajectory.
trackStates();
3652 std::vector<std::unique_ptr<GXFTrackState>> matstates;
3653 std::unique_ptr< const std::vector < const TrackStateOnSurface *>,
3654 void (*)(
const std::vector<const TrackStateOnSurface *> *) >
3656 bool matvec_used=
false;
3657 std::unique_ptr<TrackParameters> startmatpar1;
3658 std::unique_ptr<TrackParameters> startmatpar2;
3669 int npseudomuon1 = 0;
3670 int npseudomuon2 = 0;
3672 for (
auto & state : states) {
3678 if (firstidhit ==
nullptr) {
3687 if (firsthit ==
nullptr) {
3688 firsthit = state->measurement();
3690 if (tp ==
nullptr) {
3700 if (tp ==
nullptr) {
3704 state->setTrackParameters(std::unique_ptr<const TrackParameters>(tp));
3710 lasthit = state->measurement();
3716 if (firstidhit ==
nullptr) {
3717 firstidhit = state->measurement();
3720 if ((firstidpar ==
nullptr) && (tp !=
nullptr)) {
3724 lastidhit = state->measurement();
3725 if (tp !=
nullptr) {
3730 if (firstsiliconpar ==
nullptr) {
3731 firstsiliconpar = tp;
3733 lastsiliconpar = tp;
3745 if (firstmuonhit ==
nullptr) {
3746 firstmuonhit = state->measurement();
3747 if (tp !=
nullptr) {
3751 lastmuonhit = state->measurement();
3752 if (tp !=
nullptr) {
3758 if (meff->
deltaE() == 0) {
3759 if (firstcalopar ==
nullptr) {
3760 firstcalopar = state->trackParameters();
3762 lastcalopar = state->trackParameters();
3764 if (firstmatpar ==
nullptr) {
3765 firstmatpar = state->trackParameters();
3770 std::unique_ptr<TrackParameters> refpar;
3771 AmgVector(5) newpars = refpar2->parameters();
3778 newpars[0], newpars[1], newpars[2], newpars[3], newpars[4], std::nullopt
3781 if (firstmatpar !=
nullptr) {
3786 if ((startmatpar1 ==
nullptr) || ((firstidhit !=
nullptr) && (firstmuonhit !=
nullptr))) {
3790 const double mass = trajectory.
mass();
3791 if (mass > 200 * MeV) {
3792 const AmgVector(5) & newpars = startmatpar2->parameters();
3793 const double oldp = std::abs(1 / newpars[
Trk::qOverP]);
3796 startmatpar2 = startmatpar2->associatedSurface().createUniqueTrackParameters(
3797 newpars[0], newpars[1], newpars[2], newpars[3],
3798 sign / std::sqrt(oldp * oldp + 2 * 100 * MeV * std::sqrt(oldp * oldp + mass * mass) + 100 * MeV * 100 * MeV),
3803 AmgVector(5) newpars = startmatpar1->parameters();
3806 startmatpar1 = startmatpar1->associatedSurface().createUniqueTrackParameters(
3807 newpars[0], newpars[1], newpars[2], newpars[3], newpars[4], std::nullopt
3810 newpars = startmatpar2->parameters();
3813 startmatpar2 = startmatpar2->associatedSurface().createUniqueTrackParameters(
3814 newpars[0], newpars[1], newpars[2], newpars[3], newpars[4], std::nullopt
3822 refpar->momentum().unit()
3825 const double distance = getDistance(distsol);
3828 ATH_MSG_DEBUG(
"Obtaining upstream layers from Extrapolator");
3831 std::unique_ptr<const TrackParameters> tmppar;
3833 if (firstmuonhit !=
nullptr) {
3835 if (caloEntranceIsValid) {
3842 if (tmppar !=
nullptr) {
3843 destsurf = &tmppar->associatedSurface();
3848 if (matvec_used) cache.
m_matTempStore.push_back( std::move(matvec) );
3853 false, matEffects) );
3856 if (matvec && !matvec->empty()) {
3857 for (
int i = (
int)matvec->size() - 1; i > -1; i--) {
3862 std::unique_ptr<GXFMaterialEffects> meff = std::make_unique<GXFMaterialEffects>(*meot);
3863 const TrackParameters * newpars = (*matvec)[i]->trackParameters() !=
nullptr ? (*matvec)[i]->trackParameters()->clone() :
nullptr;
3864 meff->setSigmaDeltaE(0);
3865 matstates.push_back(std::make_unique<GXFTrackState>(
3867 std::unique_ptr<const TrackParameters>(newpars)
3880 refpar->momentum().unit()
3883 const double distance = getDistance(distsol);
3886 ATH_MSG_DEBUG(
"Obtaining downstream ID layers from Extrapolator");
3888 std::unique_ptr<const TrackParameters> tmppar;
3889 std::unique_ptr<Surface> calosurf;
3890 if (firstmuonhit !=
nullptr) {
3892 if (caloEntranceIsValid) {
3900 if (tmppar !=
nullptr) {
3904 cylcalosurf =
static_cast<const CylinderSurface *
>(&tmppar->associatedSurface());
3909 disccalosurf =
static_cast<const DiscSurface *
>(&tmppar->associatedSurface());
3911 if (cylcalosurf !=
nullptr) {
3914 const double radius = cylbounds.
r();
3916 calosurf = std::make_unique<CylinderSurface>(trans, radius - 1, hlength);
3917 }
else if (disccalosurf !=
nullptr) {
3918 const double newz = (
3919 disccalosurf->
center().
z() > 0 ?
3920 disccalosurf->
center().z() - 1 :
3921 disccalosurf->
center().z() + 1
3925 disccalosurf->
center().x(),
3926 disccalosurf->
center().y(),
3931 trans.translation() << newpos;
3934 const double rmin = discbounds->
rMin();
3935 const double rmax = discbounds->
rMax();
3936 calosurf = std::make_unique<DiscSurface>(trans, rmin, rmax);
3938 destsurf = calosurf.release();
3942 if (matvec_used) cache.
m_matTempStore.push_back( std::move(matvec) );
3944 ctx, *startmatpar2, *destsurf,
alongMomentum,
false, matEffects));
3945 matvec_used =
false;
3947 if (matvec && !matvec->empty()) {
3948 for (
const auto & i : *matvec) {
3954 std::unique_ptr<GXFMaterialEffects> meff = std::make_unique<GXFMaterialEffects>(*meot);
3960 meff->setDeltaE(-5);
3963 meff->setScatteringSigmas(0, 0);
3966 meff->setSigmaDeltaE(50);
3969 const TrackParameters * newparams = i->trackParameters() !=
nullptr ? i->trackParameters()->clone() :
nullptr;
3971 matstates.push_back(std::make_unique<GXFTrackState>(
3973 std::unique_ptr<const TrackParameters>(newparams)
3985 if (cache.
m_calomat && (firstmuonhit !=
nullptr) && (firstidhit !=
nullptr)) {
3988 std::vector<MaterialEffectsOnTrack> calomeots =
m_calotool->extrapolationSurfacesAndEffects(
3997 if (calomeots.empty()) {
4002 std::unique_ptr<const TrackParameters> prevtrackpars =
unique_clone(lastidpar);
4003 if (lasthit == lastmuonhit) {
4004 for (
int i = 0; i < (int) calomeots.size(); i++) {
4007 std::unique_ptr<const TrackParameters> layerpar(
m_propagator->propagateParameters(
4010 calomeots[i].associatedSurface(),
4017 if (layerpar ==
nullptr) {
4022 std::unique_ptr<GXFMaterialEffects> meff = std::make_unique<GXFMaterialEffects>(calomeots[i]);
4025 lastcalopar = layerpar.get();
4029 const double qoverp = layerpar->parameters()[
Trk::qOverP];
4030 double qoverpbrem = 0;
4034 (firstmuonpar !=
nullptr) &&
4035 std::abs(firstmuonpar->parameters()[
Trk::qOverP]) > 1.e-9
4037 qoverpbrem = firstmuonpar->parameters()[
Trk::qOverP];
4039 const double sign = (qoverp > 0) ? 1 : -1;
4040 qoverpbrem =
sign / (1 / std::abs(qoverp) - std::abs(calomeots[i].energyLoss()->deltaE()));
4043 const AmgVector(5) & newpar = layerpar->parameters();
4045 layerpar = layerpar->associatedSurface().createUniqueTrackParameters(
4046 newpar[0], newpar[1], newpar[2], newpar[3], qoverpbrem, std::nullopt
4048 meff->setdelta_p(1000 * (qoverpbrem - qoverp));
4051 matstates.push_back(std::make_unique<GXFTrackState>(
4053 std::unique_ptr<const TrackParameters>(layerpar !=
nullptr ? layerpar->clone() :
nullptr)
4055 prevtrackpars = std::move(layerpar);
4060 firsthit == firstmuonhit &&
4064 for (
int i = 0; i < (int) calomeots.size(); i++) {
4066 std::unique_ptr<const TrackParameters> layerpar(
m_propagator->propagateParameters(
4069 calomeots[i].associatedSurface(),
4076 if (layerpar ==
nullptr) {
4081 std::unique_ptr<GXFMaterialEffects> meff = std::make_unique<GXFMaterialEffects>(calomeots[i]);
4090 const double qoverpbrem = layerpar->parameters()[
Trk::qOverP];
4095 (lastmuonpar !=
nullptr) &&
4096 std::abs(lastmuonpar->parameters()[
Trk::qOverP]) > 1.e-9
4100 const double sign = (qoverpbrem > 0) ? 1 : -1;
4101 qoverp =
sign / (1 / std::abs(qoverpbrem) + std::abs(calomeots[i].energyLoss()->deltaE()));
4104 meff->setdelta_p(1000 * (qoverpbrem - qoverp));
4105 const AmgVector(5) & newpar = layerpar->parameters();
4107 prevtrackpars = layerpar->associatedSurface().createUniqueTrackParameters(
4108 newpar[0], newpar[1], newpar[2], newpar[3], qoverp, std::nullopt
4112 matstates.insert(matstates.begin(), std::make_unique<GXFTrackState>(std::move(meff), std::move(layerpar)));
4117 if (lasthit == lastmuonhit && cache.
m_extmat) {
4118 std::unique_ptr<const Trk::TrackParameters> muonpar1;
4120 if (lastcalopar !=
nullptr) {
4122 if (msEntranceIsValid) {
4130 if (muonpar1 !=
nullptr) {
4138 rot.col(2) = trackdir;
4140 trans.linear().matrix() << rot;
4141 trans.translation() << muonpar1->position() - .1 * trackdir;
4144 std::unique_ptr<const TrackParameters> curvlinpar(
m_extrapolator->extrapolateDirectly(
4152 if (curvlinpar !=
nullptr) {
4153 muonpar1 = std::move(curvlinpar);
4157 muonpar1 = std::unique_ptr<TrackParameters>(lastcalopar->
clone());
4161 muonpar1 = std::unique_ptr<TrackParameters>(refpar->clone());
4166 if (muonpar1 !=
nullptr) {
4168 muonpar1->position(),
4169 muonpar1->momentum().unit()
4173 double distance = getDistance(distsol);
4176 0) and (firstmuonhit !=
nullptr)) {
4178 muonpar1->position(),
4179 muonpar1->momentum().unit()
4185 distance = distsol.
first();
4188 std::abs(distsol.
first()) < std::abs(distsol.
second()) ?
4194 if (distance < 0 && distsol.
numberOfSolutions() > 0 && (firstidhit ==
nullptr)) {
4195 if (firstmuonpar !=
nullptr) {
4196 AmgVector(5) newpars = firstmuonpar->parameters();
4203 newpars[0], newpars[1], newpars[2], newpars[3], newpars[4], std::nullopt
4206 std::unique_ptr<const TrackParameters> tmppar(
m_propagator->propagateParameters(
4216 if (tmppar !=
nullptr) {
4217 muonpar1 = std::move(tmppar);
4223 ATH_MSG_DEBUG(
"Obtaining downstream layers from Extrapolator");
4224 if (matvec_used) cache.
m_matTempStore.push_back( std::move(matvec) );
4227 states.back()->associatedSurface(),
4231 matvec_used =
false;
4238 if (matvec && !matvec->empty()) {
4239 for (
int j = 0; j < (int) matvec->size(); j++) {
4245 std::unique_ptr<GXFMaterialEffects> meff = std::make_unique<GXFMaterialEffects>(*meot);
4250 std::abs(meff->deltaE()) > 25 &&
4251 std::abs((*matvec)[j]->trackParameters()->position().z()) < 13000
4256 const TrackParameters * newparams = (*matvec)[j]->trackParameters() !=
nullptr ? (*matvec)[j]->trackParameters()->clone() :
nullptr;
4258 matstates.push_back(std::make_unique<GXFTrackState>(
4260 std::unique_ptr<const TrackParameters>(newparams)
4270 if (firsthit == firstmuonhit && cache.
m_extmat && (firstcalopar !=
nullptr)) {
4271 std::unique_ptr<const Trk::TrackParameters> muonpar1;
4274 if (msEntranceIsValid) {
4282 if (muonpar1 !=
nullptr) {
4290 rot.col(2) = trackdir;
4292 trans.linear().matrix() << rot;
4293 trans.translation() << muonpar1->position() - .1 * trackdir;
4296 std::unique_ptr<const TrackParameters> curvlinpar(
m_extrapolator->extrapolateDirectly(
4304 if (curvlinpar !=
nullptr) {
4305 muonpar1 = std::move(curvlinpar);
4309 muonpar1 = std::unique_ptr<const TrackParameters>(firstcalopar->
clone());
4315 if (muonpar1 !=
nullptr) {
4317 muonpar1->position(),
4318 muonpar1->momentum().unit()
4322 const double distance = getDistance(distsol);
4326 ATH_MSG_DEBUG(
"Collecting upstream muon material from extrapolator");
4327 if (matvec_used) cache.
m_matTempStore.push_back( std::move(matvec) );
4330 states[0]->associatedSurface(),
4334 matvec_used =
false;
4336 if (matvec && !matvec->empty()) {
4337 ATH_MSG_DEBUG(
"Retrieved " << matvec->size() <<
" material states");
4339 for (
int j = 0; j < (int) matvec->size(); j++) {
4342 if (meb !=
nullptr) {
4347 std::unique_ptr<GXFMaterialEffects> meff = std::make_unique<GXFMaterialEffects>(*meot);
4352 std::abs(meff->deltaE()) > 25 &&
4353 std::abs((*matvec)[j]->trackParameters()->position().z()) < 13000
4359 (*matvec)[j]->trackParameters() !=
nullptr
4360 ? (*matvec)[j]->trackParameters()->clone()
4363 matstates.insert(matstates.begin(), std::make_unique<GXFTrackState>(
4365 std::unique_ptr<const TrackParameters>(tmpparams)
4378 std::vector<std::unique_ptr<GXFTrackState>> & newstates = states;
4379 std::vector<std::unique_ptr<GXFTrackState>> oldstates = std::move(newstates);
4382 newstates.reserve(oldstates.size() + matstates.size());
4385 int firstlayerno = -1;
4391 const double cosphi = std::cos(refpar->parameters()[
Trk::phi0]);
4392 const double sinphi = std::sin(refpar->parameters()[
Trk::phi0]);
4394 for (
int i = cache.
m_acceleration ? 1 : 0; i < (
int) oldstates.size(); i++) {
4395 bool addlayer =
true;
4397 while (addlayer && layerno < (
int) matstates.size()) {
4399 const TrackParameters *layerpar = matstates[layerno]->trackParameters();
4401 const DistanceSolution distsol = oldstates[i]->associatedSurface().straightLineDistanceEstimate(
4406 const double distance = getDistance(distsol);
4414 const double trackimpact = std::abs(-refpar->position().x() * sinphi + refpar->position().y() * cosphi);
4416 if (trackimpact > cylinderradius - 5 * mm) {
4422 if (i == (
int) oldstates.size() - 1) {
4431 ATH_MSG_DEBUG(
"Material state with excessive scattering, skipping it");
4442 if (firstlayerno < 0) {
4443 firstlayerno = layerno;
4448 if ((layerpar !=
nullptr) && matEffects !=
pion && matEffects !=
muon) {
4450 const Layer *lay =
nullptr;
4452 if (tvol !=
nullptr) {
4468 if (matvec_used) cache.
m_matTempStore.push_back( std::move(matvec) );
4482 const AmgVector(5) & pars = param.parameters();
4484 pars[0], pars[1], pars[2], pars[3], pars[4], std::nullopt
4493 std::unique_ptr<const TrackParameters> per(
m_extrapolator->extrapolate(
4494 Gaudi::Hive::currentContext(),param, tmppersf,
oppositeMomentum,
false, matEffects));
4496 if (per ==
nullptr) {
4497 ATH_MSG_DEBUG(
"Cannot make Perigee with starting parameters");
4501 if(std::abs(per->position().z())>5000.) {
4502 ATH_MSG_WARNING(
"Pathological perigee well outside of tracking detector!! Returning nullptr");
4510 const EventContext& ctx,
4517 ATH_MSG_DEBUG(
"--> entering GlobalChi2Fitter::myfit_helper");
4545 if (trajectory.
nDOF() < 0) {
4555 ATH_MSG_WARNING(
"Attempt to apply material corrections with q/p=0, reject track");
4569 std::unique_ptr<const TrackParameters> per =
makePerigee(cache, param, matEffects);
4585 ATH_MSG_WARNING(
"Tracking Geometry Service and/or Material Updator Tool not configured");
4596 addMaterial(ctx, cache, trajectory, per !=
nullptr ? per.get() : ¶m, matEffects);
4599 ctx, cache, trajectory, per !=
nullptr ? per.get() : ¶m, matEffects);
4611 for (
const auto & state : trajectory.
trackStates()) {
4615 state->materialEffects()->deltaE() == 0
4617 scatstate2 = state.get();
4622 scatstate = state.get();
4629 if ((scatstate ==
nullptr) || (scatstate2 ==
nullptr)) {
4630 throw std::logic_error(
"Invalid scatterer");
4635 const int nstates = (int) trajectory.
trackStates().size();
4637 trajectory.
trackStates()[nstates / 2 - 1]->position() +
4643 std::unique_ptr<const TrackParameters> nearestpar;
4644 double mindist = 99999;
4645 std::vector < GXFTrackState * >mymatvec;
4648 if ((*it).trackParameters() ==
nullptr) {
4652 const double distance = persurf
4654 (*it).trackParameters()->position(),
4655 (*it).trackParameters()->momentum().unit())
4658 const bool insideid = (
4664 (((*it).measurement() !=
nullptr) && insideid) || (
4665 ((*it).materialEffects() !=
nullptr) &&
4667 (*it).materialEffects()->deltaE() == 0 ||
4668 ((*it).materialEffects()->sigmaDeltaPhi() == 0 &&
4670 (*it).materialEffects()->deltaPhi() != 0
4674 const double dist = ((*it).trackParameters()->position() -
vertex).perp();
4675 if (dist < mindist) {
4683 if (((*it).materialEffects() !=
nullptr) && distance > 0) {
4684 mymatvec.push_back(it.get());
4688 if (nearestpar ==
nullptr) {
4692 for (
auto & state : mymatvec) {
4694 const Surface &matsurf = state->associatedSurface();
4696 nearestpar->position(), nearestpar->momentum().unit());
4698 const double distance = getDistance(distsol);
4704 std::unique_ptr<const TrackParameters> tmppar(
m_propagator->propagateParameters(
4714 if (tmppar ==
nullptr) {
4726 if (tmppar ==
nullptr) {
4736 AmgVector(5) newpars = tmppar->parameters();
4738 if (state->materialEffects()->sigmaDeltaE() > 0) {
4739 newpars[
Trk::qOverP] += .001 * state->materialEffects()->delta_p();
4742 const double de = std::abs(state->materialEffects()->deltaE());
4743 const double oldp = std::abs(1 / newpars[
Trk::qOverP]);
4744 const double newp2 = oldp * oldp - 2 * de * std::sqrt(mass * mass + oldp * oldp) + de * de;
4750 nearestpar = tmppar->associatedSurface().createUniqueTrackParameters(
4751 newpars[0], newpars[1], newpars[2], newpars[3], newpars[4], std::nullopt
4755 std::unique_ptr<Trk::TrackParameters> tmpPars(
m_propagator->propagateParameters(
4766 if (tmpPars !=
nullptr) {
4767 per.reset(
static_cast < const Perigee *
>(tmpPars.release()));
4772 const double oldp = 1. / std::abs(per->parameters()[
Trk::qOverP]);
4774 const double newp = std::sqrt(oldp * oldp + 2 * toteloss * std::sqrt(oldp * oldp + mass * mass) + toteloss * toteloss);
4775 AmgVector(5) params = per->parameters();
4778 per = per->associatedSurface().createUniqueTrackParameters(
4779 params[0], params[1], params[2], params[3], params[4], std::nullopt
4783 if (per ==
nullptr) {
4800 }
else if (per ==
nullptr) {
4817 const AmgVector(5) & pars = per->parameters();
4818 per = per->associatedSurface().createUniqueTrackParameters(
4819 pars[0], pars[1], pars[2], pars[3], 0, std::nullopt
4825 if (per !=
nullptr) {
4835 Eigen::MatrixXd a_inv;
4836 a.resize(nfitpar, nfitpar);
4841 derivPool.setZero();
4843 for (std::unique_ptr<GXFTrackState> & state : trajectory.
trackStates()) {
4844 if (state->materialEffects() !=
nullptr) {
4847 state->setDerivatives(derivPool);
4850 bool doderiv =
true;
4853 for (
int it = 0; it <
m_maxit; ++it) {
4866 runIteration(ctx, cache, trajectory, it,
a, b, lu, doderiv);
4882 const double redchi2 = (trajectory.
nDOF() > 0) ? trajectory.
chi2() / trajectory.
nDOF() : 0;
4883 const double prevredchi2 = (trajectory.
nDOF() > 0) ? trajectory.
prevchi2() / trajectory.
nDOF() : 0;
4892 (redchi2 < prevredchi2 &&
4893 (redchi2 > prevredchi2 - 1 || redchi2 < 2)) ||
4894 nsihits + ntrthits == nhits
4899 if (it != 1 || nsihits != 0 || trajectory.
nDOF() <= 0 || trajectory.
chi2() / trajectory.
nDOF() <= 3) {
4914 if (ntrtprechits+ntrttubehits) {
4915 phf = float(ntrtprechits)/float(ntrtprechits+ntrttubehits);
4917 if (phf<m_minphfcut && it>=3) {
4918 if ((ntrtprechits+ntrttubehits)>=15) {
4922 ATH_MSG_DEBUG(
"Iter = " << it <<
" | nTRTStates = " << ntrthits
4923 <<
" | nTRTPrecHits = " << ntrtprechits
4924 <<
" | nTRTTubeHits = " << ntrttubehits
4925 <<
" | nOutliers = "
4945 if (trajectory.
prefit() == 0) {
4949 Eigen::LLT < Eigen::MatrixXd >
const lltOfW(
a);
4950 if (lltOfW.info() == Eigen::Success) {
4954 const int ncols =
a.cols();
4955 Amg::MatrixX const weightInvAMG = Amg::MatrixX::Identity(ncols, ncols);
4956 a_inv = lltOfW.solve(weightInvAMG);
4975 if (traj != &trajectory) {
4980 finaltrajectory = traj;
4988 if (nperpars == 5) {
4989 for (
int i = 0; i <
a.cols(); i++) {
4990 a_inv(4, i) *= .001;
4991 a_inv(i, 4) *= .001;
4995 int scatterPos = nperpars + 2 * nscat;
4996 for (
int bremno = 0; bremno < nbrem; bremno++, scatterPos++) {
4997 for (
int i = 0; i <
a.cols(); i++) {
4998 a_inv(scatterPos, i) *= .001;
4999 a_inv(i, scatterPos) *= .001;
5006 for (
int i = 0; i < nperparams; i++) {
5007 for (
int j = 0; j < nperparams; j++) {
5008 (errmat) (j, i) = a_inv(j, i);
5013 (errmat) (4, 4) = 1e-20;
5017 std::unique_ptr<const TrackParameters> measper(
5019 perpars[0], perpars[1], perpars[2], perpars[3], perpars[4], std::move(errmat)
5029 std::unique_ptr<Track> track =
nullptr;
5031 if (finaltrajectory->
prefit() > 0) {
5032 if (finaltrajectory != &trajectory) {
5034 delete finaltrajectory;
5040 track =
makeTrack(ctx,cache, *finaltrajectory, matEffects);
5053 (track !=
nullptr) && (
5054 track->fitQuality()->numberDoF() != 0 &&
5055 track->fitQuality()->chiSquared() / track->fitQuality()->numberDoF() > cut
5058 track.reset(
nullptr);
5062 if (track ==
nullptr) {
5066 if (finaltrajectory != &trajectory) {
5067 delete finaltrajectory;
5070 return track.release();
5074 const EventContext& ctx,
5075 const Cache & cache,
5079 int & bremno_maxbrempull,
5084 std::vector<std::unique_ptr<GXFTrackState>> & states = trajectory.
trackStates();
5104 const int nmeas = (int)
res.size();
5114 const int nDOF = trajectory.
nDOF();
5115 const bool doNewPseudoMeasurements = (
5119 std::abs((trajectory.
prevchi2() - trajectory.
chi2()) / nDOF) < 15 &&
5121 (nperpars == 0 || nidhits > 0)
5132 double maxbrempull = -0.2;
5141 for (
int hitno = 0; hitno < (int) states.size(); hitno++) {
5142 std::unique_ptr<GXFTrackState> & state = states[hitno];
5156 doNewPseudoMeasurements &&
5158 !state->associatedSurface().isFree() &&
5159 !state->isRecalibrated()
5162 covMatrix(0, 0) = 100;
5164 std::unique_ptr<const PseudoMeasurementOnTrack> newpseudo = std::make_unique<const PseudoMeasurementOnTrack>(
5166 std::move(covMatrix),
5170 state->setMeasurement(std::move(newpseudo));
5171 measbase = state->measurement();
5178 double *errors = state->measurementErrors();
5180 for (
int i = 0; i < 5; i++) {
5200 res[measno] = residuals[i];
5206 res[measno] = -std::remainder(-
res[measno], 2 *
M_PI);
5215 double *errors = state->measurementErrors();
5216 for (
int i = 0; i < 5; i++) {
5217 if (errors[i] > 0) {
5218 error[measno] = errors[i];
5229 ((trajectory.
prefit() == 0) || state->materialEffects()->deltaE() == 0)
5231 const double deltaPhi = state->materialEffects()->deltaPhi();
5232 const double measDeltaPhi = state->materialEffects()->measuredDeltaPhi();
5233 const double sigma2deltaPhi = std::pow(state->materialEffects()->sigmaDeltaPhi(), 2);
5234 const double deltaTheta = state->materialEffects()->deltaTheta();
5235 const double sigma2deltaTheta = std::pow(state->materialEffects()->sigmaDeltaTheta(), 2);
5237 if (trajectory.
prefit() != 1) {
5238 b[nperpars + 2 * scatno] -= (
deltaPhi - measDeltaPhi) / sigma2deltaPhi;
5239 b[nperpars + 2 * scatno + 1] -= deltaTheta / sigma2deltaTheta;
5241 b[nperpars + scatno] -= deltaTheta / sigma2deltaTheta;
5246 deltaTheta * deltaTheta / sigma2deltaTheta
5255 if ((state->materialEffects() !=
nullptr) && state->materialEffects()->sigmaDeltaE() > 0) {
5256 double averagenergyloss = std::abs(state->materialEffects()->deltaE());
5257 const double qoverpbrem = limitInversePValue(1000 * states[hitno]->trackParameters()->parameters()[
Trk::qOverP]);
5258 const double qoverp = limitInversePValue(qoverpbrem - state->materialEffects()->delta_p());
5259 const double pbrem = 1. / std::abs(qoverpbrem);
5260 const double p = 1. / std::abs(qoverp);
5261 const double mass = .001 * trajectory.
mass();
5262 const double energy = std::sqrt(p * p + mass * mass);
5263 const double bremEnergy = std::sqrt(pbrem * pbrem + mass * mass);
5265 const double resMaterial = .001 * averagenergyloss - energy + bremEnergy;
5266 res[nmeas - nbrem + bremno] = resMaterial;
5268 const double sigde = state->materialEffects()->sigmaDeltaE();
5269 const double sigdepos = state->materialEffects()->sigmaDeltaEPos();
5270 const double sigdeneg = state->materialEffects()->sigmaDeltaENeg();
5272 double errorMaterial = .001 * state->materialEffects()->sigmaDeltaE();
5273 error[nmeas - nbrem + bremno] = errorMaterial;
5284 if (state->materialEffects()->isKink()) {
5285 maxbrempull = -999999999;
5286 state_maxbrempull =
nullptr;
5292 trajectory.
prefit() == 0 &&
5294 sigde != sigdepos &&
5297 const double elosspull = resMaterial / errorMaterial;
5299 if (trajectory.
mass() > 100) {
5304 if (std::abs(elosspull) > 1) {
5305 if (elosspull < -1) {
5306 state->materialEffects()->setSigmaDeltaE(sigdepos);
5308 state->materialEffects()->setSigmaDeltaE(sigdeneg);
5311 errorMaterial = .001 * state->materialEffects()->sigmaDeltaE();
5312 error[nmeas - nbrem + bremno] = errorMaterial;
5324 !state->materialEffects()->isKink() && (
5325 (
m_fixbrem == -1 && elosspull < maxbrempull) ||
5329 bremno_maxbrempull = bremno;
5330 state_maxbrempull = state.get();
5331 maxbrempull = elosspull;
5340 trajectory.
prefit() == 0 &&
5341 state->materialEffects()->sigmaDeltaPhi() == 0 &&
5342 state->materialEffects()->isMeasuredEloss() &&
5343 resMaterial / (.001 * state->materialEffects()->sigmaDeltaEAve()) > 2.5
5345 const TrackParameters* parforcalo = states[hitno - 2]->trackParameters();
5348 std::vector<MaterialEffectsOnTrack> calomeots =
5361 if (calomeots.size() == 3) {
5362 averagenergyloss = std::abs(calomeots[1].energyLoss()->deltaE());
5363 const double newres = .001 * averagenergyloss - energy + bremEnergy;
5364 const double newerr = .001 * calomeots[1].energyLoss()->sigmaDeltaE();
5366 const double oldPull = resMaterial / errorMaterial;
5367 const double newPull = newres / newerr;
5369 if (std::abs(newPull) < std::abs(oldPull)) {
5370 ATH_MSG_DEBUG(
"Changing from measured to parametrized energy loss");
5372 state->materialEffects()->setEloss(std::unique_ptr<EnergyLoss>(calomeots[1].energyLoss()->clone()));
5373 state->materialEffects()->setSigmaDeltaE(calomeots[1].energyLoss()->sigmaDeltaE());
5374 res[nmeas - nbrem + bremno] = newres;
5375 error[nmeas - nbrem + bremno] = newerr;
5379 state->materialEffects()->setMeasuredEloss(
false);
5389 for (
int imeas = 0; imeas < nmeas; imeas++) {
5390 if (
error[imeas] == 0) {
5405 const Cache & cache,
5411 const double oldChi2 = trajectory.
prevchi2();
5412 const double newChi2 = trajectory.
chi2();
5417 const double nDOF = trajectory.
nDOF();
5418 const double oldRedChi2 = (nDOF > 0) ? oldChi2 / nDOF : 0;
5419 const double newRedChi2 = (nDOF > 0) ? newChi2 / nDOF : 0;
5422 trajectory.
prefit() > 0 && (
5423 (newRedChi2 < 2 && it != 0) ||
5424 (newRedChi2 < oldRedChi2 + .1 && std::abs(newRedChi2 - oldRedChi2) < 1 && it != 1)
5437 int miniter = (nsihits != 0 && nsihits + ntrthits == nhits) ? 1 : 2;
5438 miniter = std::max(miniter, cache.
m_miniter);
5440 if (it >= miniter && std::abs(oldChi2 - newChi2) < 1) {
5447 const int bremno_maxbrempull,
5453 if (state_maxbrempull ==
nullptr) {
5465 const int nmeas = (int)
res.size();
5468 const double oldError =
error[nmeas - nbrem + bremno_maxbrempull];
5470 error[nmeas - nbrem + bremno_maxbrempull] = newError;
5473 if (
a.cols() != nFitPars) {
5477 const double errorRatio = oldError / newError;
5478 const double errorReductionRatio = 1 - std::pow(errorRatio, 2);
5481 for (
int i = 0; i < nFitPars; i++) {
5482 if (weightderiv(nmeas - nbrem + bremno_maxbrempull, i) == 0) {
5486 for (
int j = i; j < nFitPars; j++) {
5487 const double newaij =
a(i, j) - errorReductionRatio *
5488 weightderiv(nmeas - nbrem + bremno_maxbrempull, i) *
5489 weightderiv(nmeas - nbrem + bremno_maxbrempull, j);
5491 a.fillSymmetric(i, j, newaij);
5493 weightderiv(nmeas - nbrem + bremno_maxbrempull, i) *= errorRatio;
5502 std::vector<std::unique_ptr<GXFTrackState>> & states = trajectory.
trackStates();
5515 const int nmeas = (int) weightderiv.rows();
5517 for (std::unique_ptr<GXFTrackState> & state : states) {
5521 const auto [scatmin, scatmax] = std::minmax(scatno, nscatupstream);
5522 const auto [bremmin, bremmax] = std::minmax(bremno, nbremupstream);
5529 const double sinStereo =
5531 state->sinStereo() :
5533 const double cosStereo =
5535 std::sqrt(1 - std::pow(sinStereo, 2)) :
5543 auto getThisDeriv = [sinStereo, cosStereo, &derivatives](
int i,
int j) ->
double {
5544 if (i == 0 && sinStereo != 0) {
5545 return derivatives(0, j) * cosStereo + sinStereo * derivatives(1, j);
5547 return derivatives(i, j);
5551 for (
int i = 0; i < 5; i++) {
5566 if (i == 0 && sinStereo != 0) {
5567 weightderiv.row(measno).head(cols) =
5568 (derivatives.row(0).
head(cols) * cosStereo +
5569 sinStereo * derivatives.row(1).head(cols)) /
5572 weightderiv.row(measno).head(cols) = derivatives.row(i).head(cols) /
error[measno];
5576 for (
int j = scatmin; j < scatmax; j++) {
5577 if (trajectory.
prefit() == 1) {
5578 const int index = nperparams + j;
5581 const int index = nperparams + 2 * j;
5583 weightderiv(measno,
index + 1) = getThisDeriv(i,
index + 1) /
error[measno];
5587 for (
int j = bremmin; j < bremmax; j++) {
5588 const int index = j + nperparams + 2 * nscat;
5595 double *errors = state->measurementErrors();
5596 for (
int i = 0; i < 5; i++) {
5597 if (errors[i] > 0) {
5603 ((trajectory.
prefit() == 0) || state->materialEffects()->deltaE() == 0)
5608 if ((state->materialEffects() !=
nullptr) && state->materialEffects()->sigmaDeltaE() > 0) {
5610 const double qoverpbrem = limitInversePValue(1000 * state->trackParameters()->parameters()[
Trk::qOverP]);
5611 const double qoverp = limitInversePValue(qoverpbrem - state->materialEffects()->delta_p());
5613 const double mass = .001 * trajectory.
mass();
5615 const auto thisMeasurementIdx{nmeas - nbrem + bremno};
5619 auto multiplier = [] (
double mass,
double qOverP){
5622 const auto qoverpTerm {multiplier(mass, qoverp) /
error[thisMeasurementIdx]};
5623 const auto qoverpBremTerm {multiplier(mass, qoverpbrem) /
error[thisMeasurementIdx]};
5626 weightderiv(thisMeasurementIdx, 4) = qoverpBremTerm - qoverpTerm;
5629 const auto bremNoBase = nperparams + 2 * nscat;
5630 if (bremno < nbremupstream) {
5631 weightderiv(thisMeasurementIdx, bremNoBase + bremno) = qoverpTerm;
5632 for (
int bremno2 = bremno + 1; bremno2 < nbremupstream; bremno2++) {
5633 weightderiv(thisMeasurementIdx, bremNoBase + bremno2) = qoverpTerm - qoverpBremTerm;
5636 weightderiv(thisMeasurementIdx, bremNoBase + bremno) = qoverpBremTerm;
5637 for (
int bremno2 = nbremupstream; bremno2 < bremno; bremno2++) {
5638 weightderiv(thisMeasurementIdx, bremNoBase + bremno2) = qoverpBremTerm - qoverpTerm;
5657 const int nMeas = (int)
res.size();
5662 for (
int i = 0; i < nPerPars; i++) {
5670 for (
int i = 0; i < (int) trajectory.
trackStates().size(); i++) {
5671 const std::unique_ptr<GXFTrackState> & state = trajectory.
trackStates()[i];
5674 if (meff ==
nullptr) {
5675 measno += state->numberOfMeasuredParameters();
5679 const int firstMeasurement = i < nUpstreamStates ? 0 : measno;
5680 const int lastMeasurement = i < nUpstreamStates ? measno : nMeas - nBrem;
5683 && (trajectory.
prefit() == 0 || meff->
deltaE() == 0)) {
5684 const int scatterPos = nPerPars + 2 * scatno;
5696 const int bremPos = nPerPars + nScatPars + bremno;
5707 const Cache & cache,
5720 const int nMeas = (int)
res.size();
5722 for (
int k = 0; k < nFitPars; k++) {
5731 for (
int measno = minMeasK; measno < maxMeasK; measno++) {
5732 b[k] +=
res[measno] * (1. /
error[measno]) * weightDeriv(measno, k);
5742 if (k == 4 || k >= nPerPars + nScatPars) {
5743 for (
int measno = nMeas - nBrem; measno < nMeas; measno++) {
5744 b[k] +=
res[measno] * (1. /
error[measno]) * weightDeriv(measno, k);
5751 const Cache & cache,
5758 for (
int k = 0; k < nFitPars; k++) {
5759 for (
int l = k; l < nFitPars; l++) {
5764 for (
int measno = minMeas; measno < maxMeas; measno++) {
5765 a_kl += weightDeriv(measno, k) * weightDeriv(measno, l);
5768 a.fillSymmetric(l, k, a_kl);
5786 const int nMeas = (int)
res.size();
5793 for (
int k = nPerPars; k < nPerPars + nScatPars; k += 2) {
5794 a(k, k) += 1. / std::pow(scatSigmas[scatno].first, 2);
5795 a(k + 1, k + 1) += 1. / std::pow(scatSigmas[scatno].second, 2);
5803 for (
int measno = nMeas - nBrem; measno < nMeas; measno++) {
5804 for (
int k = 4; k < nFitPars; k++) {
5806 k = nPerPars + nScatPars;
5809 for (
int l = k; l < nFitPars; l++) {
5811 l = nPerPars + nScatPars;
5814 const double a_kl =
a(l, k) + weightDeriv(measno, k) * weightDeriv(measno, l);
5815 a.fillSymmetric(l, k, a_kl);
5827 const double oldRedChi2,
5828 const double newRedChi2
5836 bool weightChanged =
false;
5843 double newPhiWeight = 1.1;
5844 double newThetaWeight = 1.001;
5845 if (trajectory.
prefit() == 0) {
5851 newPhiWeight = 1.00000001;
5852 }
else if (it == 1) {
5853 newPhiWeight = 1.0000001;
5854 }
else if (it <= 3) {
5855 newPhiWeight = 1.0001;
5856 }
else if (it <= 6) {
5857 newPhiWeight = 1.01;
5860 if (newRedChi2 > oldRedChi2 - 1 && newRedChi2 < oldRedChi2) {
5861 newPhiWeight = 1.0001;
5862 newThetaWeight = 1.0001;
5863 }
else if (newRedChi2 > oldRedChi2 - 25 && newRedChi2 < oldRedChi2) {
5864 newPhiWeight = 1.001;
5865 newThetaWeight = 1.0001;
5872 std::size_t scatno = 0;
5877 for (
const auto & state : trajectory.
trackStates()) {
5880 if (meff ==
nullptr) {
5884 const bool isValidPlaneSurface =
5886 static_cast<const PlaneSurface *
>(&state->associatedSurface()) !=
nullptr;
5891 if (meff->
deltaE() == 0 || (trajectory.
prefit() == 0 && isValidPlaneSurface)) {
5892 weightChanged =
true;
5894 const int scatNoIndex = 2 * scatno + nPerPars;
5898 std::stringstream message;
5899 message <<
"scatno is out of range " << scatno <<
" !< " << cache.
m_phiweight.size();
5900 throw std::range_error(message.str());
5908 a(scatNoIndex, scatNoIndex) /= cache.
m_phiweight[scatno];
5912 a(scatNoIndex, scatNoIndex) *= newPhiWeight;
5913 }
else if (trajectory.
prefit() >= 2) {
5914 a(scatNoIndex, scatNoIndex) *= newPhiWeight;
5915 a(scatNoIndex + 1, scatNoIndex + 1) *= newThetaWeight;
5942 trajectory.
prefit() == 2 &&
5945 (newRedChi2 < oldRedChi2 - 25 || newRedChi2 > oldRedChi2)
5950 return weightChanged;
5959 std::size_t scatno = 0;
5969 std::stringstream message;
5970 message <<
"scatno is out of range " << scatno <<
" !< " << cache.
m_phiweight.size();
5971 throw std::range_error(message.str());
5974 const bool isValidPlaneSurface =
5976 static_cast<const PlaneSurface *
>(&state->associatedSurface()) !=
nullptr;
5978 if (meff->
deltaE() == 0 || isValidPlaneSurface) {
5979 const int scatNoIndex = 2 * scatno + nPerPars;
5980 a(scatNoIndex, scatNoIndex) /= cache.
m_phiweight[scatno];
5995 const EventContext& ctx,
6004 const int nDOFold = trajectory.
nDOF();
6005 const double oldChi2 = trajectory.
chi2();
6006 const double oldRedChi2 = nDOFold > 0 ? oldChi2 / nDOFold : 0;
6028 int bremno_maxbrempull = 0;
6044 if ((state_maxbrempull !=
nullptr) && trajectory.
converged()) {
6053 const int nDOFnew = trajectory.
nDOF();
6054 const double newChi2 = trajectory.
chi2();
6055 const double newRedChi2 = nDOFnew > 0 ? newChi2 / nDOFnew : 0;
6057 ATH_MSG_DEBUG(
"old chi2: " << oldChi2 <<
"/" << nDOFold <<
"=" << oldRedChi2 <<
6058 ", new chi2: " << newChi2 <<
"/" << nDOFnew <<
"=" << newRedChi2);
6093 if (doDeriv || weightChanged) {
6104 if (trajectory.
prefit() == 0) {
6110 if (nSiHits + nTrtHits !=
nHits) {
6117 (newRedChi2 < 2 || (newRedChi2 < oldRedChi2 && newRedChi2 > oldRedChi2 - .5))
6139 Eigen::LLT<Eigen::MatrixXd>
const llt(lu_m);
6141 if (llt.info() != Eigen::Success) {
6165 double d0 = refpar->parameters()[
Trk::d0];
6166 double z0 = refpar->parameters()[
Trk::z0];
6169 double qoverp = refpar->parameters()[
Trk::qOverP];
6171 if (nperparams > 0) {
6172 d0 += deltaParameters[0];
6173 z0 += deltaParameters[1];
6174 phi += deltaParameters[2];
6175 theta += deltaParameters[3];
6176 qoverp = (trajectory.
m_straightline) ? 0 : .001 * deltaParameters[4] + qoverp;
6188 std::vector < std::pair < double, double >>&scatangles = trajectory.
scatteringAngles();
6189 for (
int i = 0; i < nscat; i++) {
6190 scatangles[i].first += deltaParameters[2 * i + nperparams];
6191 scatangles[i].second += deltaParameters[2 * i + nperparams + 1];
6197 std::vector < double >&delta_ps = trajectory.
brems();
6198 for (
int i = 0; i < nbrem; i++) {
6199 delta_ps[i] += deltaParameters[nperparams + 2 * nscat + i];
6205 std::unique_ptr<const TrackParameters> newper(
6225 const EventContext& evtctx
6236 if (!splitProbContainer.
isValid()) {
6240 std::vector<std::unique_ptr<GXFTrackState>> & states = trajectory.
trackStates();
6247 for (
size_t stateno = 0; stateno < states.size(); stateno++) {
6252 measno += states[stateno-1]->numberOfMeasuredParameters();
6255 std::unique_ptr<GXFTrackState> & state = states[stateno];
6267 const auto *
const rot =
static_cast<const RIO_OnTrack *
>(pMeas);
6268 prd = rot->prepRawData();
6278 const auto &splitProb = splitProbContainer->splitProbability(
pixelCluster);
6279 if (!splitProb.isSplit()) {
6280 ATH_MSG_DEBUG(
"Pixel cluster is not split so no need to update" );
6284 std::unique_ptr < const RIO_OnTrack > newrot;
6285 double *olderror = state->measurementErrors();
6288 double newerror[5] = {-1,-1,-1,-1,-1};
6289 double newres[2] = {-1,-1};
6291 newrot.reset(
m_ROTcreator->correct(*prd, *trackpars, evtctx));
6296 const Amg::MatrixX & covmat = newrot->localCovariance();
6298 newerror[0] = std::sqrt(covmat(0, 0));
6299 newres[0] = newrot->localParameters()[
Trk::locX] - trackpars->parameters()[
Trk::locX];
6300 newerror[1] = std::sqrt(covmat(1, 1));
6301 newres[1] = newrot->localParameters()[
Trk::locY] - trackpars->parameters()[
Trk::locY];
6303 if (
a.cols() != nfitpars) {
6308 for(
int k =0; k<2; k++ ){
6309 const double oldres =
res[measno+k];
6310 res[measno+k] = newres[k];
6311 err[measno+k] = newerror[k];
6313 for (
int i = 0; i < nfitpars; i++) {
6314 if (weightderiv(measno+k, i) == 0) {
6318 b[i] -= weightderiv(measno+k, i) * (oldres / olderror[k] - (newres[k] * olderror[k]) / (newerror[k] * newerror[k]));
6320 for (
int j = i; j < nfitpars; j++) {
6324 weightderiv(measno+k, i) *
6325 weightderiv(measno+k, j) *
6326 ((olderror[k] * olderror[k]) / (newerror[k] * newerror[k]) - 1)
6330 weightderiv(measno+k, i) *= olderror[k] / newerror[k];
6334 state->setMeasurement(std::move(newrot));
6335 state->setMeasurementErrors(newerror);
6350 const EventContext& ctx
6358 std::vector<std::unique_ptr<GXFTrackState>> & states = trajectory.
trackStates();
6364 if (
a.cols() != nfitpars) {
6372 bool outlierremoved =
false;
6373 bool hitrecalibrated =
false;
6375 for (
int stateno = 0; stateno < (int) states.size(); stateno++) {
6376 std::unique_ptr<GXFTrackState> & state = states[stateno];
6384 std::abs(state->trackParameters()->parameters()[
Trk::driftRadius]) > 1.05 * state->associatedSurface().bounds().r()
6389 outlierremoved =
true;
6391 double *errors = state->measurementErrors();
6392 const double olderror = errors[0];
6396 for (
int i = 0; i < nfitpars; i++) {
6397 if (weightderiv(measno, i) == 0) {
6401 b[i] -=
res[measno] * weightderiv(measno, i) / olderror;
6403 for (
int j = i; j < nfitpars; j++) {
6406 a(i, j) - weightderiv(measno, i) * weightderiv(measno, j)
6409 weightderiv(measno, i) = 0;
6413 }
else if (trtrecal) {
6414 double *errors = state->measurementErrors();
6415 const double olderror = errors[0];
6417 const auto *
const thisMeasurement{state->measurement()};
6426 const double trackradius = state->trackParameters()->parameters()[
Trk::driftRadius];
6428 std::unique_ptr<const Trk::RIO_OnTrack> newrot =
nullptr;
6429 const double distance = std::abs(std::abs(trackradius) - dcradius);
6431 if (distance < scalefactor * dcerror && (olderror > 1. || trackradius * oldradius < 0)) {
6433 }
else if (distance > scalefactor * dcerror && olderror < 1.) {
6437 if (newrot !=
nullptr) {
6439 hitrecalibrated =
true;
6443 if ((measno < 0) or (measno >= (
int)
res.size())) {
6444 throw std::runtime_error(
6445 "'res' array index out of range in TrkGlobalChi2Fitter/src/GlobalChi2Fitter.cxx:" + std::to_string(__LINE__)
6449 const double oldres =
res[measno];
6450 const double newres = newradius - state->trackParameters()->parameters()[
Trk::driftRadius];
6451 errors[0] = newerror;
6452 state->setMeasurement(std::move(newrot));
6456 for (
int i = 0; i < nfitpars; i++) {
6457 if (weightderiv(measno, i) == 0) {
6461 b[i] -= weightderiv(measno, i) * (oldres / olderror - (newres * olderror) / (newerror * newerror));
6463 for (
int j = i; j < nfitpars; j++) {
6470 i < nperpars + 2 * nscats &&
6471 (i - nperpars) % 2 == 0
6478 a(i, j) + weightderiv(measno, i) * weightderiv(measno, j) * ((olderror * olderror) / (newerror * newerror) - 1) * weight
6481 weightderiv(measno, i) *= olderror / newerror;
6484 res[measno] = newres;
6485 err[measno] = newerror;
6494 measno += state->numberOfMeasuredParameters();
6498 if (trajectory.
nDOF() < 0) {
6503 if (outlierremoved || hitrecalibrated) {
6512 const EventContext& ctx,
6520 bool trackok =
false;
6522 std::unique_ptr < GXFTrajectory > cleanup_oldtrajectory;
6524 std::unique_ptr < GXFTrajectory > cleanup_newtrajectory;
6530 while (!trackok && oldtrajectory->
nDOF() > 0) {
6532 std::vector<std::unique_ptr<GXFTrackState>> & states = oldtrajectory->
trackStates();
6540 if (nhits != nsihits) {
6544 double maxsipull = -1;
6546 int hitno_maxsipull = -1;
6547 int measno_maxsipull = -1;
6548 int stateno_maxsipull = 0;
6560 for (
int stateno = 0; stateno < (int) states.size(); stateno++) {
6561 std::unique_ptr<GXFTrackState> & state = states[stateno];
6567 double *errors = state->measurementErrors();
6569 const Amg::MatrixX & hitcov = state->measurement()->localCovariance();
6570 const double sinstereo = state->sinStereo();
6571 const double cosstereo = (sinstereo == 0) ? 1 : std::sqrt(1 - sinstereo * sinstereo);
6572 double weight1 = -1;
6574 if (hitcov(0, 0) > trackcov(0, 0)) {
6575 if (sinstereo == 0) {
6576 weight1 = errors[0] * errors[0] - trackcov(0, 0);
6578 weight1 = errors[0] * errors[0] - (
6579 trackcov(0, 0) * cosstereo * cosstereo + 2 *
6580 trackcov(1, 0) * cosstereo * sinstereo + trackcov(1, 1) * sinstereo * sinstereo
6585 const double weight2 = (
6587 errors[1] * errors[1] - trackcov(1, 1) :
6591 double sipull1 = weight1 > 0 ? std::abs(
res[measno] / std::sqrt(weight1)) : -1;
6592 const double sipull2 = (
6594 std::abs(
res[measno + 1] / std::sqrt(weight2)) :
6597 sipull1 = std::max(sipull1, sipull2);
6599 if (sipull1 > maxsipull) {
6600 maxsipull = sipull1;
6601 measno_maxsipull = measno;
6602 state_maxsipull = state.get();
6603 stateno_maxsipull = stateno;
6604 hitno_maxsipull = hitno;
6615 measno += state->numberOfMeasuredParameters();
6619 const double maxpull = maxsipull;
6621 ATH_MSG_DEBUG(
" maxsipull: " << maxsipull <<
" hitno_maxsipull: " <<
6622 hitno_maxsipull <<
" n3sigma: " << n3sigma <<
" cut: " << cut <<
" cut2: " << cut2);
6633 state_maxsipull = oldtrajectory->
trackStates()[stateno_maxsipull].get();
6636 const auto *
const rot =
static_cast<const RIO_OnTrack *
>(pMeas);
6637 prd = rot->prepRawData();
6639 std::unique_ptr < const RIO_OnTrack > broadrot;
6644 Amg::VectorX parameterVector = trackpar_maxsipull->parameters();
6645 const std::unique_ptr<const TrackParameters> trackparForCorrect(
6654 state_maxsipull->trackCovariance())
6658 newerror[0] = newerror[1] = newerror[2] = newerror[3] = newerror[4] = -1;
6659 double newpull = -1;
6660 double newpull1 = -1;
6661 double newpull2 = -1;
6662 double newres1 = -1;
6663 double newres2 = -1;
6664 double newsinstereo = 0;
6677 const Amg::MatrixX & covmat = broadrot->localCovariance();
6679 if (state_maxsipull->
sinStereo() != 0) {
6680 const auto [covEigenValueSmall, covStereoAngle] = principalComponentAnalysis2x2(covmat);
6681 newerror[0] = std::sqrt(covEigenValueSmall);
6682 newsinstereo = std::sin(covStereoAngle);
6684 newerror[0] = std::sqrt(covmat(0, 0));
6687 const double cosstereo = (newsinstereo == 0) ? 1. : std::sqrt(1 - newsinstereo * newsinstereo);
6689 if (cosstereo != 1.) {
6691 cosstereo * (broadrot->localParameters()[
Trk::locX] - trackpar_maxsipull->parameters()[
Trk::locX]) +
6692 newsinstereo * (broadrot->localParameters()[
Trk::locY] - trackpar_maxsipull->parameters()[
Trk::locY])
6695 newres1 = broadrot->localParameters()[
Trk::locX] - trackpar_maxsipull->parameters()[
Trk::locX];
6698 if (newerror[0] == 0.0) {
6699 ATH_MSG_WARNING(
"Measurement error is zero or negative, treating as outlier");
6702 newpull1 = std::abs(newres1 / newerror[0]);
6706 newerror[1] = std::sqrt(covmat(1, 1));
6707 newres2 = broadrot->localParameters()[
Trk::locY] - trackpar_maxsipull->parameters()[
Trk::locY];
6708 newpull2 = std::abs(newres2 / newerror[1]);
6711 newpull = std::max(newpull1, newpull2);
6717 (newerror[0] > 1.5 * olderror[0] || newerror[1] > 1.5 * std::abs(olderror[1]))
6719 if ((measno_maxsipull < 0) or(measno_maxsipull >= (
int)
res.size())) {
6720 throw std::runtime_error(
6721 "'res' array index out of range in TrkGlobalChi2Fitter/src/GlobalChi2Fitter.cxx:" + std::to_string(__LINE__)
6726 newtrajectory = oldtrajectory;
6728 if (
a.cols() != nfitpars) {
6732 const double oldres1 =
res[measno_maxsipull];
6733 res[measno_maxsipull] = newres1;
6734 err[measno_maxsipull] = newerror[0];
6736 for (
int i = 0; i < nfitpars; i++) {
6737 if (weightderiv(measno_maxsipull, i) == 0) {
6741 b[i] -= weightderiv(measno_maxsipull, i) * (oldres1 / olderror[0] - (newres1 * olderror[0]) / (newerror[0] * newerror[0]));
6743 for (
int j = i; j < nfitpars; j++) {
6747 weightderiv(measno_maxsipull, i) *
6748 weightderiv(measno_maxsipull, j) *
6749 ((olderror[0] * olderror[0]) / (newerror[0] * newerror[0]) - 1)
6753 weightderiv(measno_maxsipull, i) *= olderror[0] / newerror[0];
6757 const double oldres2 =
res[measno_maxsipull + 1];
6758 res[measno_maxsipull + 1] = newres2;
6759 err[measno_maxsipull + 1] = newerror[1];
6761 for (
int i = 0; i < nfitpars; i++) {
6762 if (weightderiv(measno_maxsipull + 1, i) == 0) {
6766 b[i] -= weightderiv(measno_maxsipull + 1, i) * (oldres2 / olderror[1] - (newres2 * olderror[1]) / (newerror[1] * newerror[1]));
6768 for (
int j = i; j < nfitpars; j++) {
6772 weightderiv(measno_maxsipull + 1, i) *
6773 weightderiv(measno_maxsipull + 1, j) *
6774 ((olderror[1] * olderror[1]) / (newerror[1] * newerror[1]) - 1)
6779 weightderiv(measno_maxsipull + 1, i) *= olderror[1] / newerror[1];
6784 "Recovering outlier, hitno=" << hitno_maxsipull <<
" measno=" <<
6785 measno_maxsipull <<
" pull=" << maxsipull <<
" olderror_0=" <<
6786 olderror[0] <<
" newerror_0=" << newerror[0] <<
" olderror_1=" <<
6787 olderror[1] <<
" newerror_1=" << newerror[1]
6796 ((n3sigma < 2 && maxsipull > cut2 && maxsipull < cut) || n3sigma > 1) &&
6806 "Removing outlier, hitno=" << hitno_maxsipull <<
", measno=" <<
6807 measno_maxsipull <<
" pull=" << maxsipull
6814 cleanup_newtrajectory = std::make_unique<GXFTrajectory>(*oldtrajectory);
6815 newtrajectory = cleanup_newtrajectory.get();
6817 if (newa.cols() != nfitpars) {
6823 if ((measno_maxsipull < 0) or(measno_maxsipull >= (
int)
res.size())) {
6824 throw std::runtime_error(
6825 "'res' array index out of range in TrkGlobalChi2Fitter/src/GlobalChi2Fitter.cxx:" + std::to_string(__LINE__)
6829 const double oldres1 =
res[measno_maxsipull];
6830 newres[measno_maxsipull] = 0;
6832 for (
int i = 0; i < nfitpars; i++) {
6833 if (weightderiv(measno_maxsipull, i) == 0) {
6837 newb[i] -= weightderiv(measno_maxsipull, i) * oldres1 / olderror[0];
6839 for (
int j = i; j < nfitpars; j++) {
6843 weightderiv(measno_maxsipull, i) *
6844 weightderiv(measno_maxsipull, j)
6848 newweightderiv(measno_maxsipull, i) = 0;
6852 const double oldres2 =
res[measno_maxsipull + 1];
6853 newres[measno_maxsipull + 1] = 0;
6855 for (
int i = 0; i < nfitpars; i++) {
6856 if (weightderiv(measno_maxsipull + 1, i) == 0) {
6860 newb[i] -= weightderiv(measno_maxsipull + 1, i) * oldres2 / olderror[1];
6862 for (
int j = i; j < nfitpars; j++) {
6863 if (weightderiv(measno_maxsipull + 1, j) == 0) {
6870 weightderiv(measno_maxsipull + 1, i) *
6871 weightderiv(measno_maxsipull + 1, j)
6875 newweightderiv(measno_maxsipull + 1, i) = 0;
6879 newtrajectory->
setOutlier(stateno_maxsipull);
6893 for (
int it = 0; it <
m_maxit; ++it) {
6903 ctx, cache, *newtrajectory, it, *newap, *newbp, lu_m, doderiv);
6919 const double oldchi2 = oldtrajectory->
chi2() / oldtrajectory->
nDOF();
6920 const double newchi2 = (newtrajectory->
nDOF() > 0) ? newtrajectory->
chi2() / newtrajectory->
nDOF() : 0;
6923 if (newtrajectory->
nDOF() != oldtrajectory->
nDOF() && maxsipull > cut2) {
6924 mindiff = (oldchi2 > .33 *
m_chi2cut || noutl > 0) ? .8 : 1.;
6926 if (noutl == 0 && maxsipull < cut - .5 && oldchi2 < .5 *
m_chi2cut) {
6931 if (newchi2 > oldchi2 || (newchi2 > oldchi2 - mindiff && newchi2 > .33 * oldchi2)) {
6932 ATH_MSG_DEBUG(
"Outlier not confirmed, keeping old trajectory");
6940 (void)cleanup_oldtrajectory.release();
6941 return oldtrajectory;
6943 if (oldtrajectory != newtrajectory) {
6944 cleanup_oldtrajectory = std::move(cleanup_newtrajectory);
6945 oldtrajectory = newtrajectory;
6952 Eigen::LLT < Eigen::MatrixXd >
const lltOfW(
a);
6953 if (lltOfW.info() == Eigen::Success) {
6957 const int ncols =
a.cols();
6958 Amg::MatrixX const weightInvAMG = Amg::MatrixX::Identity(ncols, ncols);
6959 fullcov = lltOfW.solve(weightInvAMG);
6977 oldtrajectory->
nDOF() > 0 &&
6986 (void)cleanup_oldtrajectory.release();
6987 return oldtrajectory;
6999 if (
const auto *pMeas{hit->measurement()};
7005 nrealmeas += hit->numberOfMeasuredParameters();
7015 if (
const auto *pMeas{hit->measurement()};
7021 for (
int i = measindex; i < measindex + hit->numberOfMeasuredParameters(); i++) {
7023 cache.
m_derivmat(i, j) = derivs(measindex2, j) * errors[measindex2];
7024 if ((j == 4 && !oldtrajectory.
m_straightline) || j >= nperpars + 2 * nscat) {
7032 measindex += hit->numberOfMeasuredParameters();
7033 }
else if (hit->materialEffects() ==
nullptr) {
7034 measindex2 += hit->numberOfMeasuredParameters();
7040 const EventContext & ctx,
7048 std::unique_ptr<const TrackParameters> per(
nullptr);
7051 std::unique_ptr<const TrackParameters> prevpar(
7056 std::vector<std::pair<const Layer *, const Layer *>> & upstreamlayers = oldtrajectory.
upstreamMaterialLayers();
7059 for (
const auto & [layer1, layer2] : upstreamlayers | std::views::reverse) {
7060 if (prevpar ==
nullptr) {
7065 const Layer *layer = layer1 !=
nullptr ? layer1 : layer2;
7067 const DistanceSolution distsol = layer->surfaceRepresentation().straightLineDistanceEstimate(
7068 prevpar->position(), prevpar->momentum().unit()
7070 const double distance = getDistance(distsol);
7073 if (std::abs(distance) < 0.01) {
7077 if (distsol.
first() * distsol.
second() < 0 && !first) {
7082 if (first && distance > 0) {
7086 std::unique_ptr<const TrackParameters> layerpar(
7090 layer->surfaceRepresentation(),
7098 if (layerpar ==
nullptr) {
7102 if (layer->surfaceRepresentation().bounds().inside(layerpar->localPosition())) {
7106 prevpar = std::move(layerpar);
7119 if (startfactor > 0.5) {
7120 std::unique_ptr<const TrackParameters> updatedpar =
m_matupdator->update(
7124 if (updatedpar !=
nullptr) {
7144 if (endfactor > 0.5) {
7145 std::unique_ptr<const TrackParameters> updatedpar =
m_matupdator->update(
7149 if (updatedpar !=
nullptr) {
7155 if (prevpar !=
nullptr) {
7167 if (per ==
nullptr) {
7168 ATH_MSG_DEBUG(
"Failed to extrapolate to perigee, returning 0");
7187 std::unique_ptr<GXFTrackState>
7189 const EventContext & ctx,
7196 if (per ==
nullptr) {
7200 ATH_MSG_DEBUG(
"Final perigee: " << *per <<
" pos: " << per->position() <<
" pT: " << per->pT());
7206 const std::vector<std::unique_ptr<TrackParameters>> & hc,
7207 std::set<Identifier> & id_set,
7208 std::set<Identifier> & sct_set,
7218 for (
const std::unique_ptr<TrackParameters> & tp : hc) {
7224 if (tp ==
nullptr) {
7233 const TrkDetElementBase * de = tp->associatedSurface().associatedDetectorElement();
7235 if (de ==
nullptr) {
7246 if (id_set.find(
id) != id_set.end()) {
7265 }
else if (
m_DetID->is_sct(
id)) {
7276 }
else if (
m_DetID->is_sct(
id)) {
7285 const Identifier os = e->otherSide()->identify();
7295 if (sct_set.find(os) != sct_set.end()) {
7296 ++rv.m_sct_double_hole;
7325 for (
const std::unique_ptr<GXFTrackState> & s : trajectory.
trackStates()) {
7337 std::vector<std::reference_wrapper<GXFTrackState>> rv;
7345 for (
const std::unique_ptr<GXFTrackState> & s : trajectory.
trackStates()) {
7359 rv.emplace_back(*s);
7366 const TrkDetElementBase * de = s->trackParameters()->associatedSurface().associatedDetectorElement();
7368 if (de !=
nullptr) {
7381 if (s.get() == lastmeas) {
7391 const EventContext & ctx,
7392 const std::vector<std::reference_wrapper<GXFTrackState>> & states
7404 constexpr uint min_meas = 3;
7405 if (std::count_if(states.begin(), states.end(), [](
const GXFTrackState & s){ return s.getStateType(TrackStateOnSurface::Measurement); }) < min_meas) {
7409 bool seen_meas =
false;
7411 std::set<Identifier> id_set;
7412 std::set<Identifier> sct_set;
7418 for (std::size_t i = 0; i < states.size() - 1; i++) {
7441 const double dist = (beg.trackParameters()->position() - end.trackParameters()->position()).norm();
7443 const bool zStartValid = std::abs(beg.trackParameters()->position().z())<10000.;
7445 ATH_MSG_DEBUG(
"Pathological track parameter well outside of detector");
7446 ATH_MSG_DEBUG(
"Propagator might have issue with this, skipping");
7447 ATH_MSG_VERBOSE(
"dumping track parameters " << *(beg.trackParameters()));
7456 if (seen_meas && dist >= 2.5 && zStartValid) {
7463 std::optional<std::vector<std::unique_ptr<TrackParameters>>> & hc = beg.getHoles();
7464 std::vector<std::unique_ptr<TrackParameters>> states;
7472 if (hc.has_value()) {
7473 states = std::move(*hc);
7511 std::vector<std::unique_ptr<Trk::TrackParameters>>
const bl =
m_extrapolator->extrapolateBlindly(
7532 const EventContext & ctx,
7538 auto trajectory = std::make_unique<Trk::TrackStates>();
7546 std::unique_ptr<GXFTrackState> perigee_ts =
makeTrackFindPerigee(ctx, cache, oldtrajectory, matEffects);
7548 if (perigee_ts ==
nullptr) {
7554 trajectory->reserve(tmptrajectory.
trackStates().size());
7560 hit->resetTrackCovariance();
7566 hit->materialEffects()))
7572 auto trackState = hit->trackStateOnSurface();
7573 hit->resetTrackCovariance();
7574 trajectory->emplace_back(trackState.release());
7577 auto qual = std::make_unique<FitQuality>(tmptrajectory.
chi2(), tmptrajectory.
nDOF());
7597 std::unique_ptr<Track> rv = std::make_unique<Track>(info, std::move(trajectory), std::move(qual));
7607 std::unique_ptr<TrackSummary>
ts = std::make_unique<TrackSummary>();
7616 std::optional<TrackHoleCount> hole_count;
7623 std::vector<std::reference_wrapper<GXFTrackState>>
const states =
holeSearchStates(tmptrajectory);
7641 if (hole_count.has_value()) {
7654 rv->setTrackSummary(std::move(
ts));
7663 const EventContext & ctx,
7672 std::vector<std::unique_ptr<TrackParameters>> rv =
m_extrapolator->extrapolateStepwise(
7686 &rv.front()->associatedSurface() == &src.associatedSurface() ||
7687 trackParametersClose(*rv.front(), src, 0.001) ||
7691 rv.front().reset(
nullptr);
7702 &rv.back()->associatedSurface() == &src.associatedSurface() ||
7703 trackParametersClose(*rv.back(), src, 0.001) ||
7707 rv.back().reset(
nullptr);
7714 const EventContext & ctx,
7722 std::unique_ptr<const TrackParameters> rv;
7723 std::optional<TransportJacobian> jac{};
7734 if (rv !=
nullptr && calcderiv) {
7739 std::optional<std::vector<std::unique_ptr<TrackParameters>>> extrapolation;
7748 std::move(extrapolation)
7753 const EventContext & ctx,
7764 ctx, prev,
ts, propdir, bf, calcderiv, holesearch
7767 if (rv.m_parameters ==
nullptr) {
7768 propdir = invertPropdir(propdir);
7771 ctx, prev,
ts, propdir, bf, calcderiv, holesearch
7779 const EventContext& ctx,
7786 std::vector<std::unique_ptr<GXFTrackState>> & states = trajectory.
trackStates();
7789 std::unique_ptr<const TrackParameters> tmptrackpar;
7791 for (
int hitno = nstatesupstream - 1; hitno >= 0; hitno--) {
7792 const Surface &surf1 = states[hitno]->associatedSurface();
7799 const double distance = getDistance(distsol);
7821 (rv.m_parameters !=
nullptr) &&
7822 (prevtrackpar->
position() - rv.m_parameters->position()).mag() > 5 * mm
7828 if (rv.m_parameters ==
nullptr) {
7829 ATH_MSG_DEBUG(
"propagation failed, prev par: " << *prevtrackpar <<
7830 " pos: " << prevtrackpar->
position() <<
" destination surface: " << surf1);
7834 states[hitno]->setTrackParameters(std::move(rv.m_parameters));
7835 const TrackParameters *currenttrackpar = states[hitno]->trackParameters();
7836 const Surface &surf = states[hitno]->associatedSurface();
7838 if (rv.m_jacobian != std::nullopt) {
7840 states[hitno]->materialEffects() !=
nullptr &&
7841 states[hitno]->materialEffects()->deltaE() != 0 &&
7842 states[hitno]->materialEffects()->sigmaDeltaE() <= 0 &&
7845 const double p = 1. / std::abs(currenttrackpar->parameters()[
Trk::qOverP]);
7846 const double de = std::abs(states[hitno]->materialEffects()->deltaE());
7847 const double mass = trajectory.
mass();
7848 const double newp = std::sqrt(p * p + 2 * de * std::sqrt(mass * mass + p * p) + de * de);
7849 (*rv.m_jacobian) (4, 4) = ((p + p * de / std::sqrt(p * p + mass * mass)) / newp) * p * p / (newp * newp);
7852 states[hitno]->setJacobian(*rv.m_jacobian);
7853 }
else if (calcderiv) {
7860 if (meff !=
nullptr && hitno != 0) {
7862 surf, *meff, *states[hitno]->trackParameters(), trajectory.
mass(), -1
7865 if (std::holds_alternative<FitterStatusCode>(
r)) {
7866 return std::get<FitterStatusCode>(
r);
7869 tmptrackpar = std::move(std::get<std::unique_ptr<const TrackParameters>>(
r));
7870 prevtrackpar = tmptrackpar.get();
7872 prevtrackpar = currenttrackpar;
7878 for (
int hitno = nstatesupstream; hitno < (int) states.size(); hitno++) {
7879 const Surface &surf = states[hitno]->associatedSurface();
7883 const double distance = getDistance(distsol);
7900 (rv.m_parameters !=
nullptr) &&
7902 (prevtrackpar->
position() - rv.m_parameters->position()).mag() > 5 * mm
7907 if (rv.m_parameters ==
nullptr) {
7908 ATH_MSG_DEBUG(
"propagation failed, prev par: " << *prevtrackpar <<
7909 " pos: " << prevtrackpar->
7910 position() <<
" destination surface: " << surf);
7914 if (rv.m_jacobian != std::nullopt) {
7916 states[hitno]->materialEffects() !=
nullptr &&
7917 states[hitno]->materialEffects()->deltaE() != 0 &&
7918 states[hitno]->materialEffects()->sigmaDeltaE() <= 0 &&
7921 const double p = 1 / std::abs(rv.m_parameters->parameters()[
Trk::qOverP]);
7922 const double de = std::abs(states[hitno]->materialEffects()->deltaE());
7923 const double mass = trajectory.
mass();
7924 double newp = p * p - 2 * de * std::sqrt(mass * mass + p * p) + de * de;
7927 newp = std::sqrt(newp);
7930 (*rv.m_jacobian) (4, 4) = ((p - p * de / std::sqrt(p * p + mass * mass)) / newp) * p * p / (newp * newp);
7933 states[hitno]->setJacobian(*rv.m_jacobian);
7934 }
else if (calcderiv) {
7941 if (meff !=
nullptr) {
7943 surf, *meff, *rv.m_parameters, trajectory.
mass(), +1
7946 if (std::holds_alternative<FitterStatusCode>(
r)) {
7947 return std::get<FitterStatusCode>(
r);
7950 rv.m_parameters = std::move(std::get<std::unique_ptr<const TrackParameters>>(
r));
7953 states[hitno]->setTrackParameters(std::move(rv.m_parameters));
7954 prevtrackpar = states[hitno]->trackParameters();
7967 const AmgVector(5) & old = param.parameters();
7973 ATH_MSG_DEBUG(
"Angles out of range, phi: " << newphi <<
" theta: " << newtheta);
7977 double newqoverp = 0;
7983 const double oldp = std::abs(1 / old[
Trk::qOverP]);
7984 const double newp2 = oldp * oldp -
sign * 2 * std::abs(meff.
deltaE()) * std::sqrt(mass * mass + oldp * oldp) + meff.
deltaE() * meff.
deltaE();
7991 newqoverp = std::copysign(1 / std::sqrt(newp2), old[
Trk::qOverP]);
7998 old[0], old[1], newphi, newtheta, newqoverp, std::nullopt
8010 using Matrix55 = Eigen::Matrix<double, 5, 5>;
8012 Matrix55 initialjac;
8013 initialjac.setZero();
8014 initialjac(4, 4) = 1;
8016 Matrix55 jacvertex(initialjac);
8018 std::vector<Matrix55, Eigen::aligned_allocator<Matrix55>> jacscat(trajectory.
numberOfScatterers(), initialjac);
8019 std::vector<Matrix55, Eigen::aligned_allocator<Matrix55>> jacbrem(trajectory.
numberOfBrems(), initialjac);
8021 std::vector<std::unique_ptr<GXFTrackState>> & states = trajectory.
trackStates();
8030 for (
const bool forward : {
false,
true}) {
8032 hit_begin = nstatesupstream;
8033 hit_end = (int) states.size();
8034 scatno = nscatupstream;
8035 bremno = nbremupstream;
8037 hit_begin = nstatesupstream - 1;
8044 int hitno = hit_begin;
8045 forward ? (hitno < hit_end) : (hitno >= hit_end);
8046 hitno += (forward ? 1 : -1)
8049 state = states[hitno].get();
8053 if (fillderivmat && state->
derivatives().cols() != nfitpars) {
8061 const int jmaxbrem = 4;
8063 if (hitno == (forward ? hit_end - 1 : 0)) {
8064 if (!fillderivmat) {
8072 Eigen::Matrix<double, 5, 5> & jac = state->
jacobian();
8074 if (hitno == nstatesupstream + (forward ? 0 : -1)) {
8075 jacvertex.block<4, 5>(0, 0) = jac.block<4, 5>(0, 0);
8076 jacvertex(4, 4) = jac(4, 4);
8086 jcnt = jmax - jmin + 1;
8088 lp_bgn = forward ? nscatupstream : nscatupstream - 1;
8091 for (
int i = lp_bgn; forward ? (i < lp_end) : (i > lp_end); i += (forward ? 1 : -1)) {
8093 i == scatno + (forward ? -1 : 1) &&
8094 prevstate !=
nullptr &&
8098 jacscat[i].block(0, jmin, 4, jcnt) = jac.block(0, jmin, 4, jcnt);
8099 jacscat[i](4, 4) = jac(4, 4);
8101 calculateJac(jac, jacscat[i], jmin, jmax);
8105 Eigen::MatrixXd & derivmat = state->
derivatives();
8106 const int scatterPos = nperpars + 2 * i;
8108 derivmat.block<4, 2>(0, scatterPos) = (forward ? 1 : -1) * jacscat[i].block<4, 2>(0, 2);
8114 jcnt = jmax - jmin + 1;
8116 lp_bgn = forward ? nbremupstream : nbremupstream - 1;
8119 for (
int i = lp_bgn; forward ? (i < lp_end) : (i > lp_end); i += (forward ? 1 : -1)) {
8121 i == bremno + (forward ? -1 : 1) &&
8126 jacbrem[i].block(0, jmin, 4, jcnt) = jac.block(0, jmin, 4, jcnt);
8127 jacbrem[i](4, 4) = jac(4, 4);
8129 calculateJac(jac, jacbrem[i], jmin, jmax);
8133 Eigen::MatrixXd & derivmat = state->
derivatives();
8134 const int scatterPos = nperpars + 2 * nscats + i;
8136 derivmat.block<5, 1>(0, scatterPos) = (forward ? .001 : -.001) * jacbrem[i].block<5, 1>(0, 4);
8140 calculateJac(jac, jacvertex, 0, 4);
8144 Eigen::MatrixXd & derivmat = state->
derivatives();
8145 derivmat.block(0, 0, 4, nperpars) = jacvertex.block(0, 0, 4, nperpars);
8147 if (nperpars == 5) {
8148 derivmat.col(4).segment(0, 4) *= .001;
8149 derivmat(4, 4) = .001 * jacvertex(4, 4);
8155 (!trajectory.
prefit() || states[hitno]->materialEffects()->deltaE() == 0)
8157 scatno += (forward ? 1 : -1);
8161 states[hitno]->materialEffects() &&
8162 states[hitno]->materialEffects()->sigmaDeltaE() > 0
8164 bremno += (forward ? 1 : -1);
8167 prevstate = states[hitno].get();
8176 bool onlylocal)
const {
8182 std::vector<std::unique_ptr<GXFTrackState>> & states = trajectory.
trackStates();
8184 std::vector < int >
indices(states.size());
8186 int i = nstatesupstream;
8187 for (
int j = 0; j < (int) states.size(); j++) {
8188 if (j < nstatesupstream) {
8195 for (
int stateno = 0; stateno < (int) states.size(); stateno++) {
8196 if (stateno == 0 || stateno == nstatesupstream) {
8197 prevstate =
nullptr;
8200 std::unique_ptr<GXFTrackState> & state = states[
index];
8201 if (state->materialEffects() !=
nullptr) {
8202 prevstate = state.get();
8206 if (!state->hasTrackCovariance()) {
8207 state->zeroTrackCovariance();
8209 AmgMatrix(5, 5) & trackerrmat = state->trackCovariance();
8211 if ((prevstate !=
nullptr) &&
8215 Eigen::Matrix<double, 5, 5> & jac = state->jacobian();
8216 const AmgMatrix(5, 5)& prevcov = states[
indices[stateno - 1]]->trackCovariance();
8218 trackerrmat = jac * prevcov * jac.transpose();
8222 trackerrmat = derivatives * fullcovmat * derivatives.transpose();
8232 bool errorok =
true;
8233 for (
int i = 0; i < 5; i++) {
8236 && trackerrmat(i, i) > meascov(j, j)) {
8238 const double scale = std::sqrt(meascov(j, j) / trackerrmat(i, i));
8239 trackerrmat(i, i) = meascov(j, j);
8240 for (
int k = 0; k < 5; k++) {
8242 trackerrmat(k, i) *= scale;
8250 for (
int i = 0; i < 5; i++) {
8254 for (
int j = 0; j < 5; j++) {
8262 trackerrmat(4, 4) = 1e-20;
8266 state->trackParameters();
8268 std::optional<
AmgMatrix(5, 5)> trkerrmat;
8270 if (state->hasTrackCovariance()) {
8271 trkerrmat = (state->trackCovariance());
8273 trkerrmat = std::nullopt;
8276 const AmgVector(5) & tpars = tmptrackpar->parameters();
8277 std::unique_ptr<const TrackParameters> trackpar(
8283 std::move(trkerrmat))
8285 state->setTrackParameters(std::move(trackpar));
8288 if (errorok && trajectory.
nDOF() > 0) {
8289 fitQual =
m_updator->fullStateFitQuality(
8290 *state->trackParameters(),
8298 state->setFitQuality(fitQual);
8300 prevstate = state.get();
8304 std::optional<TransportJacobian>
8306 const EventContext& ctx,
8319 std::optional<TransportJacobian> jac = std::make_optional<TransportJacobian>(J);
8322 0.01, 0.01, 0.00001, 0.00001, 0.000000001
8333 for (
int i = 0; i < 5; i++) {
8336 if (thisdiscsurf && i == 1) {
8342 if (i == 0 && thiscylsurf) {
8343 vecminuseps[i] = -std::remainder(-vecminuseps[i], 2 *
M_PI * previousSurface.
bounds().
r());
8344 }
else if (i == 1 && thisdiscsurf) {
8345 vecpluseps[i] = -std::remainder(-vecpluseps[i], 2 *
M_PI);
8350 std::unique_ptr<const TrackParameters> parpluseps(
8360 const std::unique_ptr<const TrackParameters> parminuseps(
8371 std::unique_ptr<const TrackParameters> newparpluseps(
8382 std::unique_ptr<const TrackParameters> newparminuseps(
8397 if (newparpluseps ==
nullptr) {
8409 if (newparminuseps ==
nullptr) {
8421 if ((newparpluseps ==
nullptr) || (newparminuseps ==
nullptr)) {
8425 for (
int j = 0; j < 5; j++) {
8429 if (j == 0 && cylsurf) {
8431 }
else if (j == 1 && discsurf) {
8435 (*jac) (j, i) =
diff / (2 * eps[i]);
8448 (
"Configure the minimum number of Iterations via jobOptions");
8470 auto nmeas1 = pDataVector->size();
8471 const auto *pLastValue = (*pDataVector)[nmeas1 - 1];
8479 testrot =
static_cast<const RIO_OnTrack *
>(pLastValue);
8481 if (lastMeasIsCompetingRIO){
8483 testrot = &testcrot->rioOnTrack(0);
8487 if (testrot ==
nullptr) {
8488 const auto *pPenultimate = (*pDataVector)[nmeas1 - 2];
8491 if(penultimateIsRIO){
8492 testrot =
static_cast<const RIO_OnTrack *
>(pPenultimate);
8494 if (penultimateIsCompetingRIO){
8496 testrot = &testcrot->rioOnTrack(0);
8503 (testrot !=
nullptr) &&
8520 if (cond_obj ==
nullptr) {
8529 std::stringstream
msg;
8531 throw std::runtime_error(
msg.str());
8538 if (geometry !=
nullptr) {
8539 cache.
m_caloEntrance = geometry->trackingVolume(
"InDet::Containers::InnerDetector");
8559 if (geometry !=
nullptr) {
8560 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.