Commit 3f7e71cf authored by Johannes Josef Junggeburth's avatar Johannes Josef Junggeburth Committed by Frank Winklmeier
Browse files

STEP_Extrapolator - Insert check for valid diagonal elements

parent 1ca6b90c
......@@ -22,26 +22,39 @@
*/
namespace Amg {
/** Sometimes the extrapolation to the next surface succeeds but has termendoulsy large errors leading to uncertainties larger
than the radius of the Geneva metropole. These ones themself are clearly unphysical, but if extrapolation continues
to the next surface the numerical values blow up giving rise to floating point exception. The covariance_cutoff defines a maximum value
for the diagonal elements of the covariance matrix to consider them as sanish
*/
inline bool saneCovarianceElement(double ele) {
constexpr double upper_covariance_cutoff = 1.e34;
return !(std::isnan(ele) || std::isinf(ele) || std::abs(ele) >= upper_covariance_cutoff );
}
/** return diagonal error of the matrix
caller should ensure the matrix is symmetric and the index is in range
*/
inline double error(const Amg::MatrixX& mat, int index) {
return sqrt(mat(index, index));
return std::sqrt(mat(index, index));
}
/// Returns true if all diagonal elements of the covariance matrix
/// are greater or equal zero
template <int N> bool valid_cov(const AmgSymMatrix(N)& mat){
/// are finite, greater than zero and also below the covariance cutoff scale.
/// This check avoids floating point exceptions raised during the uncertainty calculation on the perigee parameters
/// Add a double epislon parameter to avoid numerical precision issues. The one chosen here has been
/// taken from https://stackoverflow.com/questions/1566198/how-to-portably-get-dbl-epsilon-in-c-and-c
template <int N> inline bool valid_cov(const AmgSymMatrix(N)& mat) {
static const double MIN_COV_EPSILON = 2.2204460492503131e-16;
const int dim = mat.cols();
for (int i = 0; i < dim ; ++i){
if (mat(i,i) < 0.) return false;
if ( mat(i,i) <= MIN_COV_EPSILON || !saneCovarianceElement(mat(i,i))) return false;
}
return true;
}
template<int N>
inline double error(const AmgSymMatrix(N)& mat, int index ) {
assert(index<N);
return sqrt(mat(index,index));
return std::sqrt(mat(index,index));
}
// expression template to evaluate the required size of the compressed matrix at compile time
......@@ -60,13 +73,13 @@ inline void compress(const AmgSymMatrix(N)& covMatrix, std::vector<float>& vec )
vec.reserve(CalculateCompressedSize<N>::value);
for (unsigned int i = 0; i < N ; ++i)
for (unsigned int j = 0; j <= i; ++j)
vec.push_back(covMatrix(i,j));
vec.emplace_back(covMatrix(i,j));
}
inline void compress(const MatrixX& covMatrix, std::vector<float>& vec) {
int rows = covMatrix.rows();
for (int i = 0; i < rows; ++i) {
for (int j = 0; j <= i; ++j) {
vec.push_back(covMatrix(i, j));
vec.emplace_back(covMatrix(i, j));
}
}
}
......@@ -104,10 +117,10 @@ inline void expand(std::vector<float>::const_iterator it,
template<int N>
double largestDifference(const AmgSymMatrix(N)& m1, const AmgSymMatrix(N)& m2,bool relative = false ) {
if( N < 1 ) return 0;
double max = relative ? 0.5*fabs(m1(0,0)-m2(0,0))/(m1(0,0)+m2(0,0)) : fabs(m1(0,0)-m2(0,0));
double max = relative ? 0.5*std::abs(m1(0,0)-m2(0,0))/(m1(0,0)+m2(0,0)) : std::abs(m1(0,0)-m2(0,0));
for (int i = 0; i < N; ++i) {
for (int j = 0; j < N; ++j) {
double val = fabs(m1(i,j)-m2(i,j));
double val = std::abs(m1(i,j)-m2(i,j));
if( relative ) {
val = 0.5*val/(m1(i,j)+m2(i,j));
}
......@@ -123,10 +136,10 @@ double largestDifference(const AmgSymMatrix(N)& m1, const AmgSymMatrix(N)& m2,bo
template<int N>
int largestDifference(const AmgVector(N)& m1, const AmgVector(N)& m2, bool relative = false ) {
if( N < 1 ) return 0;
double max = relative ? 0.5*fabs(m1(0)-m2(0))/(m1(0)+m2(0)) : fabs(m1(0)-m2(0));
double max = relative ? 0.5*std::abs(m1(0)-m2(0))/(m1(0)+m2(0)) : std::abs(m1(0)-m2(0));
for (int i = 0; i < N; ++i) {
for (int j = 0; j < N; ++j) {
double val = fabs(m1(i)-m2(i));
double val = std::abs(m1(i)-m2(i));
if( relative ) {
val = 0.5*val/(m1(i)+m2(i));
}
......@@ -144,9 +157,9 @@ std::pair<int, int> compare(const AmgSymMatrix(N)& m1, const AmgSymMatrix(N)& m2
for (int i = 0; i < N; ++i) {
for (int j = 0; j < N; ++j) {
if( relative ) {
if( 0.5*fabs(m1(i,j)-m2(i,j))/(m1(i,j)+m2(i,j)) > precision ) return std::make_pair(i,j);
if( 0.5*std::abs(m1(i,j)-m2(i,j))/(m1(i,j)+m2(i,j)) > precision ) return std::make_pair(i,j);
} else {
if( fabs(m1(i,j)-m2(i,j)) > precision ) return std::make_pair(i,j);
if( std::abs(m1(i,j)-m2(i,j)) > precision ) return std::make_pair(i,j);
}
}
}
......@@ -160,9 +173,9 @@ template<int N>
int compare(const AmgVector(N)& m1, const AmgVector(N)& m2, double precision = 1e-9, bool relative = false ) {
for (int i = 0; i < N; ++i) {
if( relative ) {
if( 0.5*fabs(m1(i)-m2(i))/(m1(i)+m2(i)) > precision ) return i;
if( 0.5*std::abs(m1(i)-m2(i))/(m1(i)+m2(i)) > precision ) return i;
} else {
if( fabs(m1(i)-m2(i)) > precision ) return i;
if( std::abs(m1(i)-m2(i)) > precision ) return i;
}
}
return -1;
......@@ -177,7 +190,7 @@ bool isSymMatrix(const AmgSymMatrix(N)& m ) {
if( m(i,j) < 0. ) return false;
} else {
// check that the off-diagonal elements are symmetric
if( fabs(m(i,j)-m(j,i)) > 1e-9 ) return false;
if( std::abs(m(i,j)-m(j,i)) > 1e-9 ) return false;
}
}
}
......
......@@ -147,8 +147,9 @@ namespace Muon {
bool MuonSystemExtensionTool::muonSystemExtension(const xAOD::TrackParticle& indetTrackParticle,
const MuonSystemExtension*& muonSystemExtention) const {
// get calo extension
const EventContext& ctx = Gaudi::Hive::currentContext();
std::unique_ptr<Trk::CaloExtension> caloExtension =
m_caloExtensionTool->caloExtension(Gaudi::Hive::currentContext(),
m_caloExtensionTool->caloExtension(ctx,
indetTrackParticle);
if (!caloExtension || !caloExtension->muonEntryLayerIntersection()) {
ATH_MSG_VERBOSE("Failed to get CaloExtension ");
......@@ -169,11 +170,9 @@ namespace Muon {
std::vector<std::shared_ptr<const Trk::TrackParameters> > trackParametersVec;
// loop over reference surfaces
SurfaceVec::const_iterator it = surfaces.begin();
SurfaceVec::const_iterator it_end = surfaces.end();
for (; it != it_end; ++it) {
for (const Muon::MuonLayerSurface& it : surfaces) {
// extrapolate to next layer
const Trk::Surface& surface = *it->surfacePtr;
const Trk::Surface& surface = *it.surfacePtr;
if (msgLvl(MSG::VERBOSE)) {
ATH_MSG_VERBOSE(" startPars: phi pos "
<< currentPars->position().phi() << " direction phi " << currentPars->momentum().phi() << " theta pos "
......@@ -184,33 +183,33 @@ namespace Muon {
if (currentPars->covariance())
ATH_MSG_VERBOSE(" err " << Amg::error(*currentPars->covariance(), Trk::locX) << " "
<< Amg::error(*currentPars->covariance(), Trk::locY));
ATH_MSG_VERBOSE(" destination: sector " << it->sector << " " << MuonStationIndex::regionName(it->regionIndex) << " "
<< MuonStationIndex::layerName(it->layerIndex) << " phi " << surface.center().phi()
ATH_MSG_VERBOSE(" destination: sector " << it.sector << " " << MuonStationIndex::regionName(it.regionIndex) << " "
<< MuonStationIndex::layerName(it.layerIndex) << " phi " << surface.center().phi()
<< " r " << surface.center().perp() << " z " << surface.center().z());
}
if (currentPars->momentum().perp() < 500.) {
ATH_MSG_DEBUG("Extrapolated pT less than 0.5 GeV, don't keep trying");
break;
}
const Trk::TrackParameters* exPars = m_extrapolator->extrapolate(*currentPars, surface, Trk::alongMomentum, false, Trk::muon);
std::unique_ptr<const Trk::TrackParameters> exPars {m_extrapolator->extrapolate(ctx, *currentPars, surface, Trk::alongMomentum, false, Trk::muon)};
if (!exPars) {
ATH_MSG_VERBOSE("extrapolation failed, trying next layer ");
continue;
}
// create shared pointer and add to garbage collection
std::shared_ptr<const Trk::TrackParameters> sharedPtr(exPars);
// reject intersections with very big uncertainties (parallel to surface)
if (!Amg::valid_cov(*exPars->covariance()) || Amg::error(*exPars->covariance(), Trk::locX) > 10000. ||
Amg::error(*exPars->covariance(), Trk::locY) > 10000.)
continue;
trackParametersVec.push_back(sharedPtr);
// create shared pointer and add to garbage collection
std::shared_ptr<const Trk::TrackParameters> sharedPtr{std::move(exPars)};
trackParametersVec.emplace_back(sharedPtr);
if (it->regionIndex != MuonStationIndex::Barrel) {
if (it.regionIndex != MuonStationIndex::Barrel) {
// in the endcaps we need to update the sector and check that we are in the correct sector
double phi = exPars->position().phi();
double phi = sharedPtr->position().phi();
std::vector<int> sectors;
m_sectorMapping.getSectors(phi, sectors);
......@@ -218,26 +217,26 @@ namespace Muon {
// the sector
int sector = -1;
for (auto sec : sectors) {
if (it->sector % 2 == sec % 2) {
if (it.sector % 2 == sec % 2) {
sector = sec;
break;
}
}
ATH_MSG_DEBUG(" sector " << it->sector << " selected " << sector << " nsectors " << sectors);
ATH_MSG_DEBUG(" sector " << it.sector << " selected " << sector << " nsectors " << sectors);
// only select if we found a matching sector in the set
if (sector != -1) {
MuonSystemExtension::Intersection intersection(sharedPtr, *it);
MuonSystemExtension::Intersection intersection(sharedPtr, it);
intersection.layerSurface.sector = sector;
intersections.push_back(intersection);
}
} else {
MuonSystemExtension::Intersection intersection(sharedPtr, *it);
MuonSystemExtension::Intersection intersection(sharedPtr, it);
intersections.push_back(intersection);
}
if (Amg::error(*exPars->covariance(), Trk::locX) < 10. * Amg::error(*currentPars->covariance(), Trk::locX) &&
Amg::error(*exPars->covariance(), Trk::locY) < 10. * Amg::error(*currentPars->covariance(), Trk::locY)) {
if (Amg::error(*sharedPtr->covariance(), Trk::locX) < 10. * Amg::error(*currentPars->covariance(), Trk::locX) &&
Amg::error(*sharedPtr->covariance(), Trk::locY) < 10. * Amg::error(*currentPars->covariance(), Trk::locY)) {
// only update the parameters if errors don't blow up
currentPars = exPars;
currentPars = sharedPtr.get();
}
}
......@@ -245,10 +244,7 @@ namespace Muon {
if (intersections.empty()) return false;
MuonSystemExtension* theExtension = new MuonSystemExtension(
std::unique_ptr<const Trk::TrackParameters>(caloExtension->muonEntryLayerIntersection()->clone()), std::move(intersections));
muonSystemExtention = theExtension;
muonSystemExtention = new MuonSystemExtension(caloExtension->muonEntryLayerIntersection()->uniqueClone(), std::move(intersections));
return true;
}
......
......@@ -217,7 +217,6 @@ std::unique_ptr<const Trk::TrackParameters> MuTagMatchingTool::ExtrapolateTrktoM
const Trk::TrackParameters* pTrack,
Trk::PropDirection direction) const {
if (!pSurface || !pTrack) return nullptr;
std::unique_ptr<const Trk::TrackParameters> exTrk{m_IExtrapolator->extrapolate(ctx, *pTrack, *pSurface, direction, false, Trk::muon)};
if (!exTrk) {
ATH_MSG_DEBUG(" didn't manage to extrapolate TrackParameters to abstract surface Radius " << pSurface->center().perp() << " Z "
......@@ -227,8 +226,6 @@ std::unique_ptr<const Trk::TrackParameters> MuTagMatchingTool::ExtrapolateTrktoM
}
std::unique_ptr<const Trk::Perigee> MuTagMatchingTool::flipDirection(const Trk::Perigee* inputPars) const {
// return inputPars->clone();
// CLHEP::HepVector pars = inputPars->parameters();
const AmgVector(5)& pars = inputPars->parameters();
double flippedPhi = pars[2] + M_PI;
......@@ -253,8 +250,8 @@ std::unique_ptr<const Trk::AtaPlane> MuTagMatchingTool::ExtrapolateTrktoSegmentS
Trk::PropDirection direction) const {
if (!segment || !pTrack) return nullptr;
bool isCsc(isCscSegment(segment));
unsigned int hits(cscHits(segment));
unsigned int hits{cscHits(segment)};
bool isCsc{hits>0};
if (isCsc) {
ATH_MSG_DEBUG("CSC segment has " << hits << " hits");
if (hits < 5) {
......@@ -290,7 +287,7 @@ bool MuTagMatchingTool::phiMatch(const Trk::TrackParameters* atSurface, const Mu
double PHI_CUT = m_GLOBAL_PHI_CUT;
if (hasPhi(segment)) PHI_CUT = m_GLOBAL_PHI_CUT / 2.;
const Amg::Vector3D& exTrkPos = atSurface->position();
const Amg::Vector3D exTrkPos = atSurface->position();
const Amg::Vector3D& segPos = segment->globalPosition();
const double deltaPhi = exTrkPos.deltaPhi(segPos);
......@@ -337,7 +334,7 @@ bool MuTagMatchingTool::hasPhi(const Muon::MuonSegment* seg) const {
bool MuTagMatchingTool::thetaMatch(const Trk::TrackParameters* atSurface, const Muon::MuonSegment* segment) const {
double THETA_CUT = m_GLOBAL_THETA_CUT;
const Amg::Vector3D& exTrkPos = atSurface->position();
const Amg::Vector3D exTrkPos = atSurface->position();
const Amg::Vector3D& segPos = segment->globalPosition();
const double dTheta = std::abs(exTrkPos.theta() - segPos.theta());
if (dTheta < THETA_CUT)
......@@ -349,7 +346,7 @@ bool MuTagMatchingTool::thetaMatch(const Trk::TrackParameters* atSurface, const
}
bool MuTagMatchingTool::rMatch(const Trk::TrackParameters* atSurface, const Muon::MuonSegment* segment) const {
const Amg::Vector3D& exTrkPos = atSurface->position();
const Amg::Vector3D exTrkPos = atSurface->position();
double L = exTrkPos.mag();
double R_CUT = m_GLOBAL_R_CUT * (L / 7500.); // mm
......@@ -460,7 +457,7 @@ bool MuTagMatchingTool::matchPtDependentPull(MuonCombined::MuonSegmentInfo* info
if (aMeasPer) {
double sinTheta = std::sin(aMeasPer->parameters()[Trk::theta]);
pT = sinTheta * std::abs(1.0 / (aMeasPer->parameters()[Trk::qOverP]));
pT /= 1000.;
pT *= 1.e-3;
}
double Pullcut = 2.0;
......@@ -487,8 +484,8 @@ void MuTagMatchingTool::testExtrapolation(const Trk::Surface* pSurface, const Tr
const Trk::PerigeeSurface periSurf;
std::unique_ptr<const Trk::Perigee> pPerigee =
std::make_unique<Trk::Perigee>(oripars[0], oripars[1], oripars[2], oripars[3], 0., periSurf, std::nullopt);
const Amg::Vector3D& startPos = pPerigee->position();
const Amg::Vector3D& startMom = pPerigee->momentum();
const Amg::Vector3D startPos = pPerigee->position();
const Amg::Vector3D startMom = pPerigee->momentum();
const AmgVector(5)& pars = pPerigee->parameters();
ATH_MSG_DEBUG("==============================================================================");
......@@ -502,8 +499,8 @@ void MuTagMatchingTool::testExtrapolation(const Trk::Surface* pSurface, const Tr
ATH_MSG_DEBUG("======= EXTRAPOLATED ALONG MOMENTUM ORIGINAL PERIGEE");
if (alongPars) {
const Amg::Vector3D& alongPos = alongPars->position();
const Amg::Vector3D& alongMom = alongPars->momentum();
const Amg::Vector3D alongPos = alongPars->position();
const Amg::Vector3D alongMom = alongPars->momentum();
ATH_MSG_DEBUG("=== global position " << alongPos.x() << " " << alongPos.y() << " " << alongPos.z());
ATH_MSG_DEBUG("=== global position phi theta " << alongPos.phi() << " " << alongPos.theta());
......@@ -515,28 +512,22 @@ void MuTagMatchingTool::testExtrapolation(const Trk::Surface* pSurface, const Tr
std::unique_ptr<const Trk::TrackParameters> oppositePars{
m_IExtrapolator->extrapolate(*pPerigee, *pSurface, Trk::oppositeMomentum, false, Trk::muon)};
if (oppositePars) {
const Amg::Vector3D& oppositePos = oppositePars->position();
const Amg::Vector3D& oppositeMom = oppositePars->momentum();
const Amg::Vector3D oppositePos = oppositePars->position();
const Amg::Vector3D oppositeMom = oppositePars->momentum();
ATH_MSG_DEBUG("=== global position " << oppositePos.x() << " " << oppositePos.y() << " " << oppositePos.z());
ATH_MSG_DEBUG("=== global position phi theta " << oppositePos.phi() << " " << oppositePos.theta());
ATH_MSG_DEBUG("=== global directn " << oppositeMom.phi() << " " << oppositeMom.theta());
} else
ATH_MSG_DEBUG("======= NOT EXTRAPOLATED");
double flippedPhi = pars[2] + M_PI;
if (flippedPhi > M_PI) flippedPhi -= 2 * M_PI;
double flippedTheta = M_PI - pars[3];
if (flippedTheta < 0.) flippedTheta += M_PI;
const Trk::PerigeeSurface perigSurf;
std::unique_ptr<const Trk::Perigee> flippedPerigee =
std::make_unique<Trk::Perigee>(-pars[0], pars[1], flippedPhi, flippedTheta, pars[4], perigSurf, std::nullopt);
// CLHEP::HepVector flipPars = flippedPerigee->parameters();
;
std::unique_ptr<const Trk::Perigee> flippedPerigee = flipDirection(pPerigee.get());
const AmgVector(5)& flipPars = flippedPerigee->parameters();
const Amg::Vector3D& flipPos = flippedPerigee->position();
const Amg::Vector3D& flipMom = flippedPerigee->momentum();
const Amg::Vector3D flipPos = flippedPerigee->position();
const Amg::Vector3D flipMom = flippedPerigee->momentum();
ATH_MSG_DEBUG("======= FLIPPED TRACK PARAMETERS (PERIGEE)");
ATH_MSG_DEBUG("=== phi and theta were " << pars[2] << " " << pars[3] << " and are flipped to " << flippedPhi << " " << flippedTheta);
ATH_MSG_DEBUG("=== phi and theta were " << pars[2] << " " << pars[3] << " and are flipped to ");
ATH_MSG_DEBUG("=== parameters are " << flipPars[0] << " " << flipPars[1] << " " << flipPars[2] << " " << flipPars[3] << " "
<< flipPars[4]);
ATH_MSG_DEBUG("=== global position " << flipPos.x() << " " << flipPos.y() << " " << flipPos.z());
......@@ -547,8 +538,8 @@ void MuTagMatchingTool::testExtrapolation(const Trk::Surface* pSurface, const Tr
ATH_MSG_DEBUG("======= EXTRAPOLATED ALONGFLIP MOMENTUM ORIGINAL PERIGEE");
if (alongFlipPars) {
const Amg::Vector3D& alongFlipPos = alongFlipPars->position();
const Amg::Vector3D& alongFlipMom = alongFlipPars->momentum();
const Amg::Vector3D alongFlipPos = alongFlipPars->position();
const Amg::Vector3D alongFlipMom = alongFlipPars->momentum();
ATH_MSG_DEBUG("=== global position " << alongFlipPos.x() << " " << alongFlipPos.y() << " " << alongFlipPos.z());
ATH_MSG_DEBUG("=== global position phi theta " << alongFlipPos.phi() << " " << alongFlipPos.theta());
ATH_MSG_DEBUG("=== global directn " << alongFlipMom.phi() << " " << alongFlipMom.theta());
......@@ -559,8 +550,8 @@ void MuTagMatchingTool::testExtrapolation(const Trk::Surface* pSurface, const Tr
std::unique_ptr<const Trk::TrackParameters> oppositeFlipPars{
m_IExtrapolator->extrapolate(*flippedPerigee, *pSurface, Trk::oppositeMomentum, false, Trk::muon)};
if (oppositeFlipPars) {
const Amg::Vector3D& oppositeFlipPos = oppositeFlipPars->position();
const Amg::Vector3D& oppositeFlipMom = oppositeFlipPars->momentum();
const Amg::Vector3D oppositeFlipPos = oppositeFlipPars->position();
const Amg::Vector3D oppositeFlipMom = oppositeFlipPars->momentum();
ATH_MSG_DEBUG("=== global position " << oppositeFlipPos.x() << " " << oppositeFlipPos.y() << " " << oppositeFlipPos.z());
ATH_MSG_DEBUG("=== global position phi theta " << oppositeFlipPos.phi() << " " << oppositeFlipPos.theta());
ATH_MSG_DEBUG("=== global directn " << oppositeFlipMom.phi() << " " << oppositeFlipMom.theta());
......@@ -575,12 +566,10 @@ void MuTagMatchingTool::testExtrapolation(const Trk::Surface* pSurface, const Tr
void MuTagMatchingTool::nrTriggerHits(const Muon::MuonSegment* seg, int& nRPC, int& nTGC) const {
nRPC = 0;
nTGC = 0;
std::vector<const Trk::MeasurementBase*> mbs = seg->containedMeasurements();
// for( unsigned int i = 0; i< seg->numberOfContainedROTs(); ++i){
for (unsigned int i = 0; i < mbs.size(); ++i) {
const Trk::RIO_OnTrack* rot = dynamic_cast<const Trk::RIO_OnTrack*>(mbs[i]);
for (const Trk::MeasurementBase* seg_meas : seg->containedMeasurements()) {
const Trk::RIO_OnTrack* rot = dynamic_cast<const Trk::RIO_OnTrack*>(seg_meas);
if (!rot) {
const Trk::CompetingRIOsOnTrack* crot = dynamic_cast<const Trk::CompetingRIOsOnTrack*>(mbs[i]);
const Trk::CompetingRIOsOnTrack* crot = dynamic_cast<const Trk::CompetingRIOsOnTrack*>(seg_meas);
if (crot) rot = &crot->rioOnTrack(0);
}
if (!rot) { continue; }
......@@ -828,15 +817,12 @@ MuonCombined::MuonSegmentInfo MuTagMatchingTool::muTagSegmentInfo(const Trk::Tra
info.chi2Y = chi2Y / det / 2.;
if (info.chi2Y < 0) ATH_MSG_DEBUG(" NEGATIVE chi2Y " << chi2Y << " dydyz " << dydyz << " determinant " << det);
bool hasPhi = false;
if (hitCounts.nexpectedTrigHitLayers > 1) hasPhi = true;
bool hasPhi = hitCounts.nexpectedTrigHitLayers > 1;
//
// Store hasphi
//
info.hasPhi = 0;
if (hasPhi) info.hasPhi = 1;
info.hasPhi = hasPhi;
info.t0 = 0.;
// station layer
......@@ -915,7 +901,7 @@ void MuTagMatchingTool::calculateLocalAngleErrors(const Trk::AtaPlane* exTrack,
// Parameters are described as Trk::LocX, Trk::locY, Trk::phi, Trk::theta
// So the errormatrix of the track 'localErrorMatrix' still holds global angle representation!!!!
// retrieve Jabcobian to transform the global errors err_phi,err_theta to local errors err_alphaXZ, err_alphaYZ
const Amg::RotationMatrix3D& glob2loc = exTrack->associatedSurface().transform().rotation().inverse();
const Amg::RotationMatrix3D glob2loc = exTrack->associatedSurface().transform().rotation().inverse();
const AmgVector(5)& exTrkParms = exTrack->parameters();
Trk::JacobianPhiThetaLocalAngles jacobianExTrk(exTrkParms[Trk::phi], exTrkParms[Trk::theta], glob2loc);
......@@ -966,31 +952,25 @@ bool MuTagMatchingTool::matchDistance(MuonCombined::MuonSegmentInfo* info) const
}
bool MuTagMatchingTool::isCscSegment(const Muon::MuonSegment* seg) const {
bool isCsc(false);
std::vector<const Trk::MeasurementBase*> mbs = seg->containedMeasurements();
for (unsigned int i = 0; i < mbs.size(); ++i) {
const Trk::RIO_OnTrack* rot = dynamic_cast<const Trk::RIO_OnTrack*>(mbs[i]);
for (const Trk::MeasurementBase* seg_meas : seg->containedMeasurements()) {
const Trk::RIO_OnTrack* rot = dynamic_cast<const Trk::RIO_OnTrack*>(seg_meas);
if (!rot) {
const Trk::CompetingRIOsOnTrack* crot = dynamic_cast<const Trk::CompetingRIOsOnTrack*>(mbs[i]);
const Trk::CompetingRIOsOnTrack* crot = dynamic_cast<const Trk::CompetingRIOsOnTrack*>(seg_meas);
if (crot) rot = &crot->rioOnTrack(0);
}
if (!rot) { continue; }
if (m_idHelperSvc->isCsc(rot->identify())) isCsc = true;
if (m_idHelperSvc->isCsc(rot->identify())) return true;
}
return isCsc;
return false;
}
unsigned int MuTagMatchingTool::cscHits(const Muon::MuonSegment* seg) const {
unsigned int nrHits(0);
if (!isCscSegment(seg)) return nrHits;
std::vector<const Trk::MeasurementBase*> mbs = seg->containedMeasurements();
for (unsigned int i = 0; i < mbs.size(); ++i) {
const Trk::RIO_OnTrack* rot = dynamic_cast<const Trk::RIO_OnTrack*>(mbs[i]);
unsigned int nrHits{0};
for (const Trk::MeasurementBase* seg_meas : seg->containedMeasurements()) {
const Trk::RIO_OnTrack* rot = dynamic_cast<const Trk::RIO_OnTrack*>(seg_meas);
if (!rot) {
const Trk::CompetingRIOsOnTrack* crot = dynamic_cast<const Trk::CompetingRIOsOnTrack*>(mbs[i]);
const Trk::CompetingRIOsOnTrack* crot = dynamic_cast<const Trk::CompetingRIOsOnTrack*>(seg_meas);
if (crot) rot = &crot->rioOnTrack(0);
}
if (!rot) { continue; }
......
......@@ -1385,9 +1385,15 @@ Trk::STEP_Propagator::propagateRungeKutta ( Cache&
}
//Errormatrix is included. Use Jacobian to calculate new covariance
/// Check first that the jacobian does not have crazy entries
for (int i =0;i <21;++i){
if (!Amg::saneCovarianceElement(Jacobian[i])){return nullptr;}
}
AmgSymMatrix(5) measurementCovariance = Trk::RungeKuttaUtils::newCovarianceMatrix(
Jacobian, *trackParameters->covariance());
if (!Amg::valid_cov(measurementCovariance)) return nullptr;
//Calculate multiple scattering and straggling covariance contribution.
if (cache.m_matPropOK && (m_multipleScattering || m_straggling) && std::abs(totalPath)>0.) {
if (returnCurv || targetSurfaces[solutions[0]].first->type()==Trk::SurfaceType::Cone) {
......
Markdown is supported
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment