Muon::MuonTrackCleaner Class Reference

#include <MuonTrackCleaner.h>

Inheritance diagram for Muon::MuonTrackCleaner:

Collaboration diagram for Muon::MuonTrackCleaner:

Classes
struct	ChamberLayerStatistics
struct	ChamberPullInfo
struct	ChamberRemovalOutput
	struct to store return values of chamber removal, contains the new track plus a list the removed hits More...
struct	CleaningState
struct	EtaPhiHits
struct	MCTBCleaningInfo
struct	SortByAvePull
struct	SortChamberRemovalResultByChi2Ndof

Public Types
typedef std::vector< MCTBCleaningInfo >	InfoVec
typedef InfoVec::iterator	InfoIt
typedef std::set< const Trk::MeasurementBase * >	MeasSet
typedef MeasSet::iterator	MeasIt
typedef std::map< Identifier, ChamberPullInfo >	PullChamberMap
typedef PullChamberMap::iterator	PullChamberIt
typedef PullChamberMap::const_iterator	PullChamberCit
typedef std::map< Identifier, EtaPhiHits >	EtaPhiPerChamberMap
typedef EtaPhiPerChamberMap::iterator	EtaPhiPerChamberIt
typedef std::vector< std::pair< double, Identifier > >	PullChVec
typedef PullChVec::iterator	PullChIt

Public Member Functions
	~MuonTrackCleaner ()=default
	destructor More...
StatusCode	initialize () override
	AlgTool initialize. More...
std::unique_ptr< Trk::Track >	clean (const Trk::Track &track, const EventContext &ctx) const override
	clean a track, returns a pointer to a new track if successfull. More...
std::unique_ptr< Trk::Track >	clean (const Trk::Track &track, const std::set< Identifier > &chamberRemovalExclusionList, const EventContext &ctx) const override
	clean a track, returns a pointer to a new track if successfull. More...

Private Member Functions
const Trk::ResidualPull *	calculateResPul (const Trk::MeasurementBase &meas, const Trk::TrackParameters &pars) const
	calculate Residual/Pull for a given MeasurementBase + TrackParameters, ownership is transfered to user More...
std::unique_ptr< Trk::Track >	cleanTrack (const EventContext &ctx, const Trk::Track *track, CleaningState &state) const
	clean a track, actual implementation More...
double	calcPull (const double residual, const double locMesCov, const double locTrkCov, const bool &trkStateIsUnbiased) const
	calculate the pull given measurement error and track error More...
std::unique_ptr< Trk::Track >	cleanCompROTs (const EventContext &ctx, std::unique_ptr< Trk::Track > track, CleaningState &state) const
	clean up competing ROTs that consist out of two clusters More...
std::unique_ptr< Trk::Track >	recoverFlippedMdt (const EventContext &ctx, std::unique_ptr< Trk::Track > track, CleaningState &state) const
	flip signs of MDT hits with large pull if pull if the oppositely signed radius is small More...
std::unique_ptr< Trk::Track >	hitCleaning (const EventContext &ctx, std::unique_ptr< Trk::Track > track, CleaningState &state) const
	remove bad hits from track. More...
std::unique_ptr< Trk::Track >	chamberCleaning (const EventContext &ctx, std::unique_ptr< Trk::Track > track, CleaningState &state) const
	remove bad chamber from track. More...
std::unique_ptr< Trk::Track >	outlierRecovery (const EventContext &ctx, std::unique_ptr< Trk::Track > track, CleaningState &state, const MuonStationIndex::ChIndex *currentIndex=nullptr) const
	recover outliers that are within the cuts. More...
bool	isOutsideOnTrackCut (const Identifier &id, double res, double pull, double cutScaleFactor) const
	check whether hit is an outlier More...
ChamberRemovalOutput	removeChamber (const EventContext &ctx, const std::unique_ptr< Trk::Track > &track, Identifier chId, bool removePhi, bool removeEta, CleaningState &state) const
	remove chamber from track More...
void	init (const EventContext &ctx, const Trk::Track &track, CleaningState &state) const
	init cleaner More...
bool	checkStations (CleaningState &state) const
bool	checkInnerConstraint (CleaningState &state) const
bool	checkPhiConstraint (CleaningState &state) const
void	printStates (Trk::Track *track) const
std::unique_ptr< Trk::Track >	fitTrack (const EventContext &ctx, Trk::Track &track, Trk::ParticleHypothesis pHyp, bool slFit) const
bool	extractChambersToBeRemoved (CleaningState &state, std::set< Identifier > &chambersToBeRemovedSet, bool usePhi=false) const
	helper function to extract chambers that are to be removed More...

Static Private Member Functions
static void	unremoveHits (ChamberRemovalOutput &result)
static std::string	print (ChamberLayerStatistics &statistics)

Private Attributes
ToolHandle< Trk::ITrackFitter >	m_trackFitter
ToolHandle< Trk::ITrackFitter >	m_slTrackFitter {this, "SLFitter", "Trk::GlobalChi2Fitter/MCTBSLFitterMaterialFromTrack"}
ToolHandle< Trk::IUpdator >	m_measurementUpdator {this, "MeasurementUpdator", "Trk::KalmanUpdator/MuonMeasUpdator"}
ToolHandle< Muon::IMdtDriftCircleOnTrackCreator >	m_mdtRotCreator
ToolHandle< IMuonCompetingClustersOnTrackCreator >	m_compRotCreator
ToolHandle< Trk::IResidualPullCalculator >	m_pullCalculator
ServiceHandle< Muon::IMuonEDMHelperSvc >	m_edmHelperSvc
PublicToolHandle< Muon::MuonEDMPrinterTool >	m_printer {this, "Printer", "Muon::MuonEDMPrinterTool/MuonEDMPrinterTool"}
ServiceHandle< Muon::IMuonIdHelperSvc >	m_idHelperSvc {this, "MuonIdHelperSvc", "Muon::MuonIdHelperSvc/MuonIdHelperSvc"}
ToolHandle< Trk::IExtrapolator >	m_extrapolator {this, "Extrapolator", "Trk::Extrapolator/AtlasExtrapolator"}
SG::ReadCondHandleKey< AtlasFieldCacheCondObj >	m_fieldCacheCondObjInputKey
Gaudi::Property< bool >	m_useMdtResiCut {this, "UseMdtResiCut", false}
Gaudi::Property< double >	m_chi2Cut {this, "Chi2Cut", 100.}
Gaudi::Property< double >	m_pullCut {this, "PullCut", 5.}
Gaudi::Property< double >	m_mdtResiCut {this, "MdtResiCut", 1.}
Gaudi::Property< double >	m_pullCutPhi {this, "PullCutPhi", 10.}
Gaudi::Property< double >	m_associationScaleFactor {this, "AssociationScaleFactor", 0.7}
Gaudi::Property< unsigned int >	m_ncycles {this, "CleaningCycles", 5}
Gaudi::Property< double >	m_avePullSumPerChamberCut {this, "MaxAvePullSumPerChamber", 3.5}
Gaudi::Property< bool >	m_recoverOutliers {this, "RecoverOutliers", true}
Gaudi::Property< bool >	m_flipMdtDriftRadii {this, "FlipMdtDriftRadii", true}
Gaudi::Property< bool >	m_cleanCompROTs {this, "CleanCompROTs", true}
Gaudi::Property< bool >	m_onlyUseHitErrorInRecovery {this, "OnlyUseHitErrorInRecovery", true}
Gaudi::Property< double >	m_adcCut {this, "AdcCut", 50.}
Gaudi::Property< bool >	m_iterate {this, "Iterate", 0.7}
Gaudi::Property< bool >	m_use_slFit {this, "UseSLFit", false}

Detailed Description

Definition at line 44 of file MuonTrackCleaner.h.

Member Typedef Documentation

EtaPhiPerChamberIt

typedef EtaPhiPerChamberMap::iterator Muon::MuonTrackCleaner::EtaPhiPerChamberIt

Definition at line 129 of file MuonTrackCleaner.h.

EtaPhiPerChamberMap

typedef std::map<Identifier, EtaPhiHits> Muon::MuonTrackCleaner::EtaPhiPerChamberMap

Definition at line 128 of file MuonTrackCleaner.h.

InfoIt

typedef InfoVec::iterator Muon::MuonTrackCleaner::InfoIt

Definition at line 94 of file MuonTrackCleaner.h.

InfoVec

typedef std::vector<MCTBCleaningInfo> Muon::MuonTrackCleaner::InfoVec

Definition at line 93 of file MuonTrackCleaner.h.

MeasIt

typedef MeasSet::iterator Muon::MuonTrackCleaner::MeasIt

Definition at line 97 of file MuonTrackCleaner.h.

MeasSet

typedef std::set<const Trk::MeasurementBase*> Muon::MuonTrackCleaner::MeasSet

Definition at line 96 of file MuonTrackCleaner.h.

PullChamberCit

typedef PullChamberMap::const_iterator Muon::MuonTrackCleaner::PullChamberCit

Definition at line 121 of file MuonTrackCleaner.h.

PullChamberIt

typedef PullChamberMap::iterator Muon::MuonTrackCleaner::PullChamberIt

Definition at line 120 of file MuonTrackCleaner.h.

PullChamberMap

typedef std::map<Identifier, ChamberPullInfo> Muon::MuonTrackCleaner::PullChamberMap

Definition at line 119 of file MuonTrackCleaner.h.

PullChIt

typedef PullChVec::iterator Muon::MuonTrackCleaner::PullChIt

Definition at line 132 of file MuonTrackCleaner.h.

PullChVec

typedef std::vector<std::pair<double, Identifier> > Muon::MuonTrackCleaner::PullChVec

Definition at line 131 of file MuonTrackCleaner.h.

Constructor & Destructor Documentation

~MuonTrackCleaner()

Muon::MuonTrackCleaner::~MuonTrackCleaner ( )

default

destructor

Member Function Documentation

calcPull()

double Muon::MuonTrackCleaner::calcPull ( const double  residual, const double  locMesCov, const double  locTrkCov, const bool &  trkStateIsUnbiased  ) const

private

calculate the pull given measurement error and track error

calculateResPul()

const Trk::ResidualPull* Muon::MuonTrackCleaner::calculateResPul ( const Trk::MeasurementBase &  meas, const Trk::TrackParameters &  pars  ) const

private

calculate Residual/Pull for a given MeasurementBase + TrackParameters, ownership is transfered to user

chamberCleaning()

std::unique_ptr< Trk::Track > Muon::MuonTrackCleaner::chamberCleaning ( const EventContext &  ctx, std::unique_ptr< Trk::Track >  track, CleaningState &  state  ) const

private

remove bad chamber from track.

The returned track pointer can be zero, pointing to the initial track or a new pointer

Definition at line 462 of file MuonTrackCleaner.cxx.

                                                                                              {
        ATH_MSG_DEBUG("run chamber cleaning on track " << m_printer->print(*track));
 
        if (state.chambersToBeRemoved.empty() && state.chambersToBeRemovedPhi.empty()) return track;
 
        if (state.pullSumPerChamber.size() == 2) {
            // loop over chambers that should be removed, if one of them is MDT or CSC don't clean
            // as this would result in a single chamber track
            PullChIt chit = state.chambersToBeRemoved.begin();
            PullChIt chit_end = state.chambersToBeRemoved.end();
            for (; chit != chit_end; ++chit) {
                const Identifier& chid = chit->second;
                if (m_idHelperSvc->isMdt(chid) || m_idHelperSvc->isCsc(chid)) {
                    ATH_MSG_DEBUG(" only two precision chambers, cannot remove chamber.  ");
                    return nullptr;
                }
            }
        }
 
        if (!state.chamberRemovalExclusionList.empty()) {
            unsigned int foundChambers = 0;
            // loop over chambers that should be removed and check if they are in the exclusion list
            PullChIt chit = state.chambersToBeRemoved.begin();
            PullChIt chit_end = state.chambersToBeRemoved.end();
            for (; chit != chit_end; ++chit) {
                const Identifier& chid = chit->second;
                if (m_idHelperSvc->isMdt(chid) && state.chamberRemovalExclusionList.count(chid)) {
                    ATH_MSG_DEBUG(" found excluded chamber " << m_idHelperSvc->toStringChamber(chid));
                    ++foundChambers;
                }
            }
            unsigned int excludedChambers = state.chamberRemovalExclusionList.size();
            if (foundChambers > excludedChambers) {
                ATH_MSG_WARNING(" Found more excluded chambers than in list, this should not happen ");
                return nullptr;
            } else if (foundChambers == excludedChambers) {
                ATH_MSG_DEBUG(" all excluded chambers in removal list, failing cleaning ");
                return nullptr;
            }
        }
 
        std::vector<ChamberRemovalOutput> cleaningResults;
        cleaningResults.reserve(state.chambersToBeRemoved.size());
 
        // TODO first remove hits in both eta and phi
 
        ATH_MSG_DEBUG(" number of chambers short listed to be removed: " << state.chambersToBeRemoved.size());
        std::stable_sort(state.chambersToBeRemoved.begin(), state.chambersToBeRemoved.end(), SortByAvePull());
 
        // try all combinations of removing a chamber with eta hits
        unsigned int nchambers = state.chambersToBeRemoved.size() + state.chambersToBeRemovedPhi.size();
        PullChIt chit = state.chambersToBeRemoved.begin();
        PullChIt chit_end = state.chambersToBeRemoved.end();
        for (; chit != chit_end; ++chit) {
            ChamberRemovalOutput result = removeChamber(ctx, track, chit->second, false, true, state);
            if (!result.track) {
                ATH_MSG_DEBUG(" Removed eta hits of " << m_idHelperSvc->toStringChamber(chit->second) << ", track lost ");
                continue;
            }
            ATH_MSG_DEBUG(" Removed eta hits of " << m_idHelperSvc->toStringChamber(chit->second));
            ATH_MSG_DEBUG(m_printer->print(*result.track));
 
            // this is an optimization for the common case where there is only one chamber to be removed
            if (nchambers > 1) unremoveHits(result);
 
            cleaningResults.push_back(std::move(result));
        }  // for (chit)
 
        ATH_MSG_DEBUG(" number of phi chambers short listed to be removed: " << state.chambersToBeRemoved.size());
        std::stable_sort(state.chambersToBeRemovedPhi.begin(), state.chambersToBeRemovedPhi.end(), SortByAvePull());
        // try all combinations of removing a chamber with phi hits
        chit = state.chambersToBeRemovedPhi.begin();
        chit_end = state.chambersToBeRemovedPhi.end();
        for (; chit != chit_end; ++chit) {
            ChamberRemovalOutput result = removeChamber(ctx, track, chit->second, true, false, state);
            if (!result.track) {
                ATH_MSG_DEBUG(" Removed phi hits of " << m_idHelperSvc->toStringChamber(chit->second) << ", track lost ");
                continue;
            }
            ATH_MSG_DEBUG(" Removed phi hits of " << m_idHelperSvc->toStringChamber(chit->second));
            ATH_MSG_DEBUG(m_printer->print(*result.track));
 
            // this is an optimization for the common case where there is only one chamber to be removed
            if (nchambers > 1) unremoveHits(result);
 
            cleaningResults.push_back(std::move(result));
        }
 
        // if no cleaned tracks return 0
        if (cleaningResults.empty()) return nullptr;
        std::stable_sort(cleaningResults.begin(), cleaningResults.end(), SortChamberRemovalResultByChi2Ndof());
        ATH_MSG_DEBUG(" chamberCleaning Results nr " << cleaningResults.size());
        for (auto& res : cleaningResults) { ATH_MSG_DEBUG(" track " << m_printer->print(*res.track)); }
 
        ChamberRemovalOutput& finalResult = cleaningResults.front();
        if (nchambers > 1) {
            // loop over removed hits and redo 'remove'
            for (auto* hit : finalResult.removedHits) hit->useInFit = 0;
        }
 
        ATH_MSG_DEBUG(" chamberCleaning:  track ");
        ATH_MSG_DEBUG(m_printer->print(*finalResult.track));
 
        // make clone just in case outlier recovery fails
        std::unique_ptr<Trk::Track> finalResultTrackClone = std::make_unique<Trk::Track>(*finalResult.track);
 
        init(ctx, *finalResultTrackClone, state);
 
        // now finally check whether the removed layer now is recoverable (happens sometimes if the segment has one or more bad hits)
        ChIndex removedChamberIndex = m_idHelperSvc->chamberIndex(finalResult.chId);
        std::unique_ptr<Trk::Track> recoveredTrack = outlierRecovery(ctx, std::move(finalResult.track), state, &removedChamberIndex);
        if (!recoveredTrack) return finalResultTrackClone;
        init(ctx, *recoveredTrack, state);
        return recoveredTrack;
    }

checkInnerConstraint()

bool Muon::MuonTrackCleaner::checkInnerConstraint ( CleaningState &  state ) const

private

Definition at line 1519 of file MuonTrackCleaner.cxx.

                                                                          {
        unsigned int nstations = state.stations.size();
        if (nstations == 1 ||
            (nstations == 2 && (state.stations.count(StIndex::EM) && state.stations.count(StIndex::EO)))) {
            ATH_MSG_DEBUG(" Momentum measurement lost, cleaning given up ");
            return false;
        }
        return true;
    }

checkPhiConstraint()

bool Muon::MuonTrackCleaner::checkPhiConstraint ( CleaningState &  state ) const

private

Definition at line 1529 of file MuonTrackCleaner.cxx.

                                                                        {
        if (state.stations.size() == 1 && !state.nIdHits && state.nPhiConstraints < 2) {
            ATH_MSG_DEBUG(" Underconstraint fit " << state.nPhiConstraints);
            return false;
        }
        return true;
    }

checkStations()

bool Muon::MuonTrackCleaner::checkStations ( CleaningState &  state ) const

private

Definition at line 1508 of file MuonTrackCleaner.cxx.

                                                                   {
        unsigned int nstationsChamberCleaning = state.stations.size();
        ATH_MSG_DEBUG(" Cleaner removed full station from track: remaining layers " << nstationsChamberCleaning);
 
        if (nstationsChamberCleaning < 2) {
            ATH_MSG_DEBUG(" Cleaner removed all but one station from track!!! ");
            return false;
        }
        return true;
    }

clean() [1/2]

std::unique_ptr< Trk::Track > Muon::MuonTrackCleaner::clean ( const Trk::Track &  track, const EventContext &  ctx  ) const

override

clean a track, returns a pointer to a new track if successfull.

If the input track is does not require cleaning a pointer the the initial track is return in which case the user should not delete the old track! The caller should ensure the track gets deleted.

Definition at line 65 of file MuonTrackCleaner.cxx.

                                                                                                        {
        CleaningState state;
        std::unique_ptr<Trk::Track> cleanedTrack = cleanTrack(ctx, &track, state);
        return cleanedTrack;
    }

clean() [2/2]

std::unique_ptr< Trk::Track > Muon::MuonTrackCleaner::clean ( const Trk::Track &  track, const std::set< Identifier > &  chamberRemovalExclusionList, const EventContext &  ctx  ) const

override

clean a track, returns a pointer to a new track if successfull.

If the input track is does not require cleaning a pointer the the initial track is return in which case the user should not delete the old track! The cleaning will not clean if all the chambers in the exclusions list are marked as to be deleted. The caller should ensure the track gets deleted.

Definition at line 53 of file MuonTrackCleaner.cxx.

                                                                                       {
        CleaningState state;
        state.chamberRemovalExclusionList = chamberRemovalExclusionList;
 
        if (!state.chamberRemovalExclusionList.empty()) {
            ATH_MSG_DEBUG(" Cleaning with exclusion list " << state.chamberRemovalExclusionList.size());
        }
 
        std::unique_ptr<Trk::Track> cleanedTrack = cleanTrack(ctx, &track, state);
        return cleanedTrack;
    }

cleanCompROTs()

std::unique_ptr< Trk::Track > Muon::MuonTrackCleaner::cleanCompROTs ( const EventContext &  ctx, std::unique_ptr< Trk::Track >  track, CleaningState &  state  ) const

private

clean up competing ROTs that consist out of two clusters

Definition at line 174 of file MuonTrackCleaner.cxx.

                                                                                            {
        if (!m_cleanCompROTs || state.numberOfCleanedCompROTs == 0) return track;
 
        const Trk::Perigee* perigee = track->perigeeParameters();
        if (!perigee) {
            ATH_MSG_DEBUG("   track without perigee ");
            return nullptr;
        }
 
        ATH_MSG_DEBUG(" Clean comp rots " << state.numberOfCleanedCompROTs);
 
        auto tsos = std::make_unique<Trk::TrackStates>();
        tsos->reserve(state.measInfo.size());
 
        unsigned int nmeas = 0;
        // loop over hits
        InfoIt hit = state.measInfo.begin();
        InfoIt hit_end = state.measInfo.end();
        for (; hit != hit_end; ++hit) {
            // hits that are flagged as outlier or hits in the chamber to be removed are added as Outlier
            if (!hit->useInFit) {
                ATH_MSG_DEBUG("   removing hit " << m_idHelperSvc->toString(hit->id) << " pull " << hit->resPull->pull().front());
                if (hit->inBounds)
                    tsos->push_back(MuonTSOSHelper::cloneTSOSWithUpdate(*hit->originalState, *hit->meas, *hit->pars,
                                                                       Trk::TrackStateOnSurface::Outlier));
 
                continue;
            } else {
                if (hit->meas) ++nmeas;
 
                if (hit->cleanedCompROT) {
                    tsos->push_back(MuonTSOSHelper::cloneTSOSWithUpdate(*hit->originalState, *hit->cleanedCompROT, *hit->pars,
                                                                       Trk::TrackStateOnSurface::Measurement));
                    ATH_MSG_DEBUG("   replacing CompROT " << m_idHelperSvc->toString(hit->id) << " pull " << hit->resPull->pull().front());
                } else {
                    tsos->push_back(hit->originalState->clone());
                }
            }
        }
 
        if (nmeas < 6) {
            ATH_MSG_DEBUG(" too few hits, cannot recover CompROTS ");
            return nullptr;
        }
 
        // create new track
        std::unique_ptr<Trk::Track> cleanedTrack =
            std::make_unique<Trk::Track>(track->info(), std::move(tsos), track->fitQuality() ? track->fitQuality()->uniqueClone() : nullptr);
        printStates(cleanedTrack.get());
 
        // fit new track
        std::unique_ptr<Trk::Track> newTrack = fitTrack(ctx, *cleanedTrack, track->info().particleHypothesis(), state.slFit);
 
        if (newTrack) {
            init(ctx, *newTrack, state);
            return newTrack;
        } else
            return nullptr;
    }

cleanTrack()

std::unique_ptr< Trk::Track > Muon::MuonTrackCleaner::cleanTrack ( const EventContext &  ctx, const Trk::Track *  track, CleaningState &  state  ) const

private

clean a track, actual implementation

Definition at line 71 of file MuonTrackCleaner.cxx.

                                                                                                                                   {
        ATH_MSG_DEBUG(" Perform cleaning for track ");
        ATH_MSG_DEBUG(m_printer->print(*track));
 
        init(ctx, *track, state);
 
        ATH_MSG_DEBUG("after init, track is " << m_printer->print(*track));
        ATH_MSG_DEBUG("  start cleaning ");
 
        unsigned int nstationsInitial = state.stations.size();
 
        // first clean up chambers
        std::unique_ptr<Trk::Track> chamberTrack = chamberCleaning(ctx, std::make_unique<Trk::Track>(*track), state);
        if (!chamberTrack) {
            ATH_MSG_DEBUG(" chamber removal failed ");
            return nullptr;
        }
        ATH_MSG_DEBUG("after chamber cleaning, track is " << m_printer->print(*chamberTrack));
 
        // fail if only station remains on the track
        unsigned int nstationsChamberCleaning = state.stations.size();
        if (nstationsInitial != nstationsChamberCleaning) {
            if (!checkStations(state)) return nullptr;
        }
 
        // if performed curved fit and we did not have a vertex constraint or ID hits, reject the track if inner was removed
        if (!state.slFit && !state.hasVertexConstraint && state.nIdHits <= 0) {
            if (!checkInnerConstraint(state)) return nullptr;
        }
 
        // if performing a single station layer cleaning without ID hits, reject track if there are insufficient phi constraints
        if (!checkPhiConstraint(state)) return nullptr;
 
        // clean competing ROTs
        std::unique_ptr<Trk::Track> cleanCompTrack = cleanCompROTs(ctx, std::move(chamberTrack), state);
        if (!cleanCompTrack) {
            ATH_MSG_DEBUG(" CompROT cleaning failed ");
            return nullptr;
        }
        ATH_MSG_DEBUG("after comp rot cleaning, track is " << m_printer->print(*cleanCompTrack));
 
        // recover MDTs with flipped signs
        std::unique_ptr<Trk::Track> flippedTrack = recoverFlippedMdt(ctx, std::move(cleanCompTrack), state);
        if (!flippedTrack) {
            ATH_MSG_DEBUG(" MDT sign flipping failed ");
            return nullptr;
        }
        ATH_MSG_DEBUG("after flipped mdt recovery, track is " << m_printer->print(*flippedTrack));
 
        std::unique_ptr<Trk::Track> hitTrack = hitCleaning(ctx, std::move(flippedTrack), state);
        if (!hitTrack) {
            ATH_MSG_DEBUG(" track lost after outlier removal ");
            return nullptr;
        }
        ATH_MSG_DEBUG("after hit cleaning, track is " << m_printer->print(*hitTrack));
 
        // keep this clone in case outlier recovery fails but the rest of the cleaning was succesful.
        std::unique_ptr<Trk::Track> hitTrackClone = std::make_unique<Trk::Track>(*hitTrack);
 
        // fail if only one station remains
        unsigned int nstationsHitCleaning = state.stations.size();
        if (nstationsInitial != nstationsHitCleaning) {
            if (!checkStations(state)) return nullptr;
        }
 
        // if performed curved fit and we did not have a vertex constraint or ID hits, reject the track if inner was removed
        if (!state.slFit && !state.hasVertexConstraint && state.nIdHits <= 0) {
            if (!checkInnerConstraint(state)) return nullptr;
        }
 
        // if performing a single station layer cleaning without ID hits, reject track if there are insufficient phi constraints
        if (!checkPhiConstraint(state)) return nullptr;
 
        std::unique_ptr<Trk::Track> cleanedTrack = outlierRecovery(ctx, std::move(hitTrack), state);
        // do not discard tracks that fail outlierRecovery, check that the track is ok
        // note that this also performs a useful check on the quality of the cleaning in general
        if (!cleanedTrack || !state.chambersToBeRemoved.empty() || !state.largePullMeasurements.empty()) {
            init(ctx, *hitTrackClone, state);
            if (!state.chambersToBeRemoved.empty() || !state.largePullMeasurements.empty()) {
                ATH_MSG_DEBUG("Outlier recovery failure unrecoverable, reject track");
                return nullptr;
            } else {
                ATH_MSG_DEBUG("Outlier recovery failed but initial track is recoverable");
                cleanedTrack = std::move(hitTrackClone);
            }
        }
 
        if (state.nhits < state.noutliers) {
            ATH_MSG_DEBUG(" track rejected due to high outlier ratio: hits " << state.nhits << "   outliers " << state.noutliers);
            return nullptr;
        }
 
        unsigned int nstationsFinal = state.stations.size();
        if (nstationsInitial != nstationsFinal) {
            if (!checkStations(state)) return nullptr;
        }
 
        ATH_MSG_VERBOSE(" final track " << m_printer->print(*cleanedTrack));
        ATH_MSG_VERBOSE(m_printer->printStations(*cleanedTrack));
 
        return cleanedTrack;
    }

extractChambersToBeRemoved()

bool Muon::MuonTrackCleaner::extractChambersToBeRemoved ( CleaningState &  state, std::set< Identifier > &  chambersToBeRemovedSet, bool  usePhi = false  ) const

private

helper function to extract chambers that are to be removed

Definition at line 1402 of file MuonTrackCleaner.cxx.

                                                                         {
        bool doCleaning = false;
        PullChamberMap& pullSumPerChamber = usePhi ? state.pullSumPerChamberPhi : state.pullSumPerChamber;
        PullChVec& chambersToBeRemoved = usePhi ? state.chambersToBeRemovedPhi : state.chambersToBeRemoved;
 
        PullChamberCit cit = pullSumPerChamber.begin();
        PullChamberCit cit_end = pullSumPerChamber.end();
        if (msgLvl(MSG::DEBUG)) {
            if (!pullSumPerChamber.empty()) msg() << MSG::DEBUG << "Chamber pulls " << pullSumPerChamber.size() << ":";
        }
        int ndof = 0;
        double pulltot = 0.;
        for (; cit != cit_end; ++cit) {
            double avePull = cit->second.pullSum / cit->second.nhits;
            pulltot += cit->second.pullSum;
            ndof += cit->second.nhits;
            double avePullReduced = cit->second.nhits > 1 ? (cit->second.pullSum - cit->second.maxPull) / (cit->second.nhits - 1) : avePull;
            if (msgLvl(MSG::DEBUG))
                msg() << MSG::DEBUG << std::endl
                      << " chamber " << m_idHelperSvc->toStringChamber(cit->first) << "  pull sum " << cit->second.pullSum << " max pull "
                      << cit->second.maxPull << " nhits " << cit->second.nhits << " ave pull " << avePull << " reduced ave pull "
                      << avePullReduced;
            if (avePull > m_avePullSumPerChamberCut) {
                doCleaning = true;
                if (cit->second.nhits > 3 && avePullReduced < m_avePullSumPerChamberCut) {
                    if (msgLvl(MSG::DEBUG)) msg() << MSG::DEBUG << "  large pull sum ignored ";
                    doCleaning = false;
                } else {
                    if (msgLvl(MSG::DEBUG)) msg() << MSG::DEBUG << "  large pull sum => removing";
                    if (!chambersToBeRemovedSet.count(cit->first)) {
                        chambersToBeRemoved.emplace_back(avePull, cit->first);
                        chambersToBeRemovedSet.insert(cit->first);
                    }
                }
            }
        }
        if (msgLvl(MSG::DEBUG)) msg() << endmsg;
        bool dropMore = false;
        if (dropMore && pulltot * pulltot > 2. * ndof * ndof) {
            doCleaning = true;
            if (msgLvl(MSG::DEBUG)) msg() << MSG::DEBUG << "  large pull per dof " << pulltot << " ndof " << ndof << endmsg;
        }
 
        if (doCleaning && m_iterate) {
            cit = pullSumPerChamber.begin();
            for (; cit != cit_end; ++cit) {
                double avePull = cit->second.pullSum / cit->second.nhits;
                if (!chambersToBeRemovedSet.count(cit->first)) {
                    chambersToBeRemoved.emplace_back(avePull, cit->first);
                    chambersToBeRemovedSet.insert(cit->first);
                }
            }
        }
 
        if (dropMore && doCleaning && m_iterate) {
            // add trigger chambers
            cit = state.pullSumPerChamberEta.begin();
            cit_end = state.pullSumPerChamberEta.end();
            for (; cit != cit_end; ++cit) {
                double avePull = cit->second.pullSum / cit->second.nhits;
                if (!chambersToBeRemovedSet.count(cit->first)) {
                    chambersToBeRemoved.emplace_back(avePull, cit->first);
                    chambersToBeRemovedSet.insert(cit->first);
                }
            }
        }
 
        return doCleaning;
    }

fitTrack()

std::unique_ptr< Trk::Track > Muon::MuonTrackCleaner::fitTrack ( const EventContext &  ctx, Trk::Track &  track, Trk::ParticleHypothesis  pHyp, bool  slFit  ) const

private

Definition at line 1549 of file MuonTrackCleaner.cxx.

                                                                             {
        return slFit ? m_slTrackFitter->fit(ctx, track, false, pHyp) : m_trackFitter->fit(ctx, track, false, pHyp);
    }

hitCleaning()

std::unique_ptr< Trk::Track > Muon::MuonTrackCleaner::hitCleaning ( const EventContext &  ctx, std::unique_ptr< Trk::Track >  track, CleaningState &  state  ) const

private

remove bad hits from track.

The returned track pointer can be zero, pointing to the initial track or a new pointer

Definition at line 296 of file MuonTrackCleaner.cxx.

                                                                                          {
        if (state.largePullMeasurements.empty()) return track;
        ATH_MSG_DEBUG(" trying outlier removal ");
 
        const Trk::Perigee* perigee = track->perigeeParameters();
        if (!perigee) {
            ATH_MSG_DEBUG("input track without perigee ");
            return nullptr;
        }
 
        std::unique_ptr<Trk::Track> newTrack;
 
        for (unsigned int n = 0; n < m_ncycles; ++n) {
            // sanity check, should not remove too many hits
            if (state.largePullMeasurements.size() > 10) {
                ATH_MSG_DEBUG(" Too many outliers, cannot perform cleaning ");
                return nullptr;
            }
 
            ATH_MSG_VERBOSE(" outlier removal cycle " << n);
 
            auto tsos = std::make_unique<Trk::TrackStates>();
            tsos->reserve(state.measInfo.size());
            ATH_MSG_VERBOSE("cleaning track with " << state.measInfo.size() << " hits");
 
            unsigned int nmeas = 0;
            unsigned int nremovedPhi = 0;
            bool hasSmall = false;
            bool hasLarge = false;
            MCTBCleaningInfo* firstPhi = nullptr;
            MCTBCleaningInfo* lastPhi = nullptr;
            std::map<StIndex, std::pair<bool, bool> > slCountsPerStationLayer;
            // loop over hits
            InfoIt hit = state.measInfo.begin();
            InfoIt hit_end = state.measInfo.end();
            for (; hit != hit_end; ++hit) {
                bool remove = state.largePullMeasurements.count(hit->meas) || !hit->inBounds || hit->isNoise;
                // hits that are flagged as outlier or hits in the chamber to be removed are added as Outlier
                if (!hit->useInFit || remove) {
                    hit->useInFit = 0;
 
                    // count number of phi outliers
                    if (!hit->id.is_valid() || m_idHelperSvc->measuresPhi(hit->id)) ++nremovedPhi;
 
                    if (remove) {
                        std::string inb = hit->inBounds ? " inBounds" : " outBounds";
                        ATH_MSG_DEBUG("   removing hit " << m_idHelperSvc->toString(hit->id) << " pull " << hit->resPull->pull().front()
                                                         << inb);
                    }
                    if (hit->inBounds) {
                        if (hit->cleanedCompROT) {
                            tsos->push_back(MuonTSOSHelper::cloneTSOSWithUpdate(*hit->originalState, *hit->cleanedCompROT, *hit->pars,
                                                                               Trk::TrackStateOnSurface::Outlier));
                        } else {
                            tsos->push_back(MuonTSOSHelper::cloneTSOSWithUpdate(*hit->originalState, *hit->meas, *hit->pars,
                                                                               Trk::TrackStateOnSurface::Outlier));
                        }
                    }
                    continue;
                } else {
                    if (hit->resPull)
                        ATH_MSG_DEBUG("   keeping hit " << m_idHelperSvc->toString(hit->id) << " pull " << hit->resPull->pull().front());
 
                    if (hit->meas) {
                        ++nmeas;
                        if (!hit->id.is_valid() || m_idHelperSvc->measuresPhi(hit->id)) {
                            if (!firstPhi)
                                firstPhi = &(*hit);
                            else
                                lastPhi = &(*hit);
                        }
 
                        if (hit->id.is_valid() && m_idHelperSvc->isMdt(hit->id)) {
                            StIndex stIndex = toStationIndex(hit->chIndex);
                            bool isSmall = m_idHelperSvc->isSmallChamber(hit->id);
                            bool isLarge = !isSmall;
                            // look for layer
                            std::map<StIndex, std::pair<bool, bool> >::iterator pos =
                                slCountsPerStationLayer.find(stIndex);
 
                            // if not found add current
                            if (pos == slCountsPerStationLayer.end())
                                slCountsPerStationLayer[stIndex] = std::make_pair(isSmall, isLarge);
                            else {
                                // update
                                if (isSmall) slCountsPerStationLayer[stIndex].first = true;
                                if (isLarge) slCountsPerStationLayer[stIndex].second = true;
                            }
                        }
                    }
                    if (!hit->originalState) ATH_MSG_DEBUG("no original state!");
                    tsos->push_back(hit->originalState->clone());
                }
            }
            // loop over sl map and count the overlaps
            unsigned int noverlaps = 0;
            for (auto [index, slCounts] : slCountsPerStationLayer) {
                if (!hasSmall && slCounts.first) hasSmall = true;
                if (!hasLarge && slCounts.second) hasLarge = true;
                if (slCounts.first && slCounts.second) ++noverlaps;
            }
            // now update for the case there are small and large chambers but not in the same layer
            if (noverlaps == 0 && hasSmall && hasLarge) ++noverlaps;
 
            if (nmeas < 6) {
                ATH_MSG_DEBUG(" too few hits, cannot perform hitCleaning ");
                return nullptr;
            }
 
            ATH_MSG_DEBUG(" nremovedPhi " << nremovedPhi << " noverlaps " << noverlaps << " nid " << state.nIdHits);
            if (firstPhi) ATH_MSG_DEBUG(" hasFirstPhi: " << m_idHelperSvc->toString(firstPhi->id));
            if (lastPhi) ATH_MSG_DEBUG(" hasLastPhi: " << m_idHelperSvc->toString(lastPhi->id));
 
            // only perform check on phi constraints if any phi hits were removed
            if (nremovedPhi > 0) {
                bool hasPhiConstraint = false;
 
                if (state.nIdHits > 0)
                    hasPhiConstraint = true;  // ok if ID hits
                else if (firstPhi && noverlaps > 0)
                    hasPhiConstraint = true;  // ok if one phi hit + one overlap
                else if (noverlaps > 1)
                    hasPhiConstraint = true;  // ok if two overlaps
                else if (firstPhi && lastPhi && firstPhi->pars && lastPhi->pars) {
                    double distPhi =
                        std::abs((firstPhi->pars->position() - lastPhi->pars->position()).dot(firstPhi->pars->momentum().unit()));
                    ATH_MSG_DEBUG(" Distance between phi hits " << distPhi);
                    if (distPhi > 450.) hasPhiConstraint = true;
                }
                if (!hasPhiConstraint) {
                    ATH_MSG_DEBUG("Lost phi constraint during track cleaning, reject track");
                    return nullptr;
                }
            }
 
            // create new track
            std::unique_ptr<Trk::Track> cleanedTrack =
                std::make_unique<Trk::Track>(track->info(), std::move(tsos), track->fitQuality() ? track->fitQuality()->uniqueClone() : nullptr);
 
            // fit new track
            printStates(cleanedTrack.get());
 
            newTrack = fitTrack(ctx, *cleanedTrack, track->info().particleHypothesis(), state.slFit);
 
            if (!newTrack) {
                return nullptr;
            } else {
                const Trk::Perigee* perigee = newTrack->perigeeParameters();
                if (!perigee) {
                    ATH_MSG_DEBUG("   track without perigee ");
                    return nullptr;
                }
                // reinitialize cleaner
                init(ctx, *newTrack, state);
            }
 
            if (state.largePullMeasurements.empty()) {
                ATH_MSG_DEBUG("   cleaning ended successfully after cycle " << n);
                return newTrack;
            }
        }
 
        return nullptr;
    }

init()

void Muon::MuonTrackCleaner::init ( const EventContext &  ctx, const Trk::Track &  track, CleaningState &  state  ) const

private

init cleaner

Definition at line 788 of file MuonTrackCleaner.cxx.

                                                                                                          {
        // Tell clang to optimize assuming that FP exceptions can trap.
        // Otherwise, it can vectorize the division, which can lead to
        // spurious division-by-zero traps from unused vector lanes.
        CXXUTILS_TRAPPING_FP;
 
        state.nscatterers = 0;
        state.numberOfFlippedMdts = 0;
        state.numberOfCleanedCompROTs = 0;
        state.nhits = 0;
        state.noutliers = 0;
        state.pullSum = ChamberPullInfo();
        state.pullSumPhi = ChamberPullInfo();
        state.pullSumTrigEta = ChamberPullInfo();
        state.chambersToBeRemoved.clear();
        state.chambersToBeRemovedPhi.clear();
        state.pullSumPerChamber.clear();
        state.pullSumPerChamberPhi.clear();
        state.pullSumPerChamberEta.clear();
        state.hitsPerChamber.clear();
        state.outBoundsPerChamber.clear();
        state.measInfo.clear();
        state.largePullMeasurements.clear();
        state.chamberLayerStatistics.clear();
        state.stations.clear();
        state.phiLayers.clear();
        state.hasOfBoundsOutliers = false;
        state.hasVertexConstraint = false;
        state.hasSmall = false;
        state.hasLarge = false;
        state.nIdHits = 0;
        state.nPseudoMeasurements = 0;
        state.nPhiConstraints = 0;
 
        MagField::AtlasFieldCache fieldCache;
        // Get field cache object
        SG::ReadCondHandle<AtlasFieldCacheCondObj> readHandle{m_fieldCacheCondObjInputKey, ctx};
        if (!readHandle.isValid()) {
            ATH_MSG_ERROR("Failed to retrieve AtlasFieldCacheCondObj with key " << m_fieldCacheCondObjInputKey.key());
            return;
        }
        readHandle->getInitializedCache(fieldCache);
 
        state.slFit = !fieldCache.toroidOn() || m_edmHelperSvc->isSLTrack(track);
        if (m_use_slFit) state.slFit = true;
 
        // loop over track and calculate residuals
        const Trk::TrackStates* states = track.trackStateOnSurfaces();
        if (!states) {
            ATH_MSG_WARNING(" track without states, cannot perform cleaning ");
            return;
        }
 
        ATH_MSG_DEBUG(" init: " << m_printer->print(track));
 
        state.measInfo.reserve(states->size());
 
        std::set<int> rpcLayers;
        std::set<int> tgcLayers;
        const Trk::MeasurementBase* mdtmeas = nullptr;
        double largestmdtpull = -999;
 
        // loop over TSOSs
        for (const Trk::TrackStateOnSurface* tsit : *states) {
            if (tsit->type(Trk::TrackStateOnSurface::Perigee))
                ATH_MSG_DEBUG(tsit->dumpType() << ", parameters " << *tsit->trackParameters());
            const Trk::TrackParameters* pars = tsit->trackParameters();
            if (!pars) {
                state.measInfo.emplace_back(tsit);
                continue;
            }
 
            if (tsit->type(Trk::TrackStateOnSurface::Scatterer)) {
                const Trk::MaterialEffectsBase* matEff = tsit->materialEffectsOnTrack();
                const Trk::MaterialEffectsOnTrack* matTrk = dynamic_cast<const Trk::MaterialEffectsOnTrack*>(matEff);
                ATH_MSG_VERBOSE(" Scatterer: r  " << pars->position().perp() << " z " << pars->position().z());
                if (matEff) ATH_MSG_VERBOSE(" X0 " << matEff->thicknessInX0());
                if (matTrk && matTrk->scatteringAngles()) {
                    const Trk::ScatteringAngles* scatAngle = matTrk->scatteringAngles();
                    ATH_MSG_VERBOSE(" pull phi " << scatAngle->deltaPhi() / scatAngle->sigmaDeltaPhi() << " pull theta "
                                                 << scatAngle->deltaTheta() / scatAngle->sigmaDeltaTheta());
                }
                if (matTrk && matTrk->energyLoss()) {
                    const Trk::EnergyLoss* eloss = matTrk->energyLoss();
                    ATH_MSG_DEBUG(" eloss " << eloss->deltaE());
                }
                ++state.nscatterers;
                state.measInfo.emplace_back(tsit);
                continue;
            }
 
            // check whether state is a measurement
            const Trk::MeasurementBase* meas = tsit->measurementOnTrack();
            if (!meas) {
                state.measInfo.emplace_back(tsit);
                continue;
            }
 
            Identifier id = m_edmHelperSvc->getIdentifier(*meas);
            bool pseudo = !id.is_valid();
            if (pseudo) ++state.nPseudoMeasurements;
 
            if (!pseudo && !m_idHelperSvc->mdtIdHelper().is_muon(id)) {
                ATH_MSG_VERBOSE("      TSOS is not a muon hit, position: r  " << pars->position().perp() << " z " << pars->position().z());
 
                // count ID hits on track
                if (tsit->type(Trk::TrackStateOnSurface::Measurement) && m_idHelperSvc->mdtIdHelper().is_indet(id)) { ++state.nIdHits; }
                state.measInfo.emplace_back(tsit);
                continue;
            }
 
            // check whether there is a vertex constraint
            if (pseudo && dynamic_cast<const Trk::PseudoMeasurementOnTrack*>(meas) && pars->associatedSurface().center().perp() < 200.) {
                state.hasVertexConstraint = true;
            }
 
            Identifier chId = pseudo ? id : m_idHelperSvc->chamberId(id);
            bool measuresPhi = pseudo ? true : m_idHelperSvc->measuresPhi(id);
 
            // bound checks
            Amg::Vector2D locPos;
            if (!meas->associatedSurface().globalToLocal(pars->position(), pars->position(), locPos)) {
                ATH_MSG_DEBUG(" localToGlobal failed !!!!! ");
                continue;
            }
            bool inBounds = true;
            double tol1 = 100.;
            double tol2 = 2 * tol1;
            if (!pseudo && m_idHelperSvc->isMdt(id)) tol1 = 5.;
 
            // we need a special bound check for MDTs so we cast to SL surface
            const Trk::StraightLineSurface* slSurf = dynamic_cast<const Trk::StraightLineSurface*>(&meas->associatedSurface());
            // we need a special bound check also for MMs to consider edge passivation
            const MMClusterOnTrack* mmClusterOnTrack = dynamic_cast<const MMClusterOnTrack*>(meas);
 
            if (slSurf) {
                // perform bound check only for second coordinate
                inBounds = slSurf->bounds().insideLoc2(locPos, tol2);
            } else if (mmClusterOnTrack) {
                // for MM, perform the bound check from the detector element
                inBounds = mmClusterOnTrack->detectorElement()->insideActiveBounds(id, locPos, tol1, tol2);
            } else {
                inBounds = meas->associatedSurface().insideBounds(locPos, tol1, tol2);
            }
 
            ChIndex chIndex = !pseudo ? m_idHelperSvc->chamberIndex(id) : ChIndex::ChUnknown;
 
            // pointer to resPull: workaround because a const pointer is returned
            std::optional<Trk::ResidualPull> resPull{m_pullCalculator->residualPull(
                meas, pars, tsit->type(Trk::TrackStateOnSurface::Outlier) ? Trk::ResidualPull::Unbiased : Trk::ResidualPull::Biased)};
            if (!resPull) {
                ATH_MSG_DEBUG(" calculation of residual/pull failed !!!!! ");
                continue;
            }
            int pullSize = resPull->pull().size();
            double residual = resPull->residual().front();
            double pull = std::abs(resPull->pull().front());
 
            // sanity check
            if (!pseudo && pullSize != 1) {
                ATH_MSG_WARNING(" ResidualPull vector has size " << pullSize << " for channel " << m_idHelperSvc->toString(id));
                continue;
            }
 
            bool isMDT = !pseudo ? m_idHelperSvc->isMdt(id) : false;
            double error = pull > 0.001 ? std::abs(residual / pull) : 1000.;
            double rDrift = isMDT ? meas->localParameters()[Trk::locR] : 0.;
            double rTrack = isMDT ? pars->parameters()[Trk::locR] : 0.;
            double rTrackAbs = std::abs(rTrack);
            double flippedResidual = isMDT ? rDrift + rTrack : 1e10;
            double flippedPull = isMDT ? std::abs(flippedResidual / error) : 1e10;
 
            bool isNoise = false;
            bool isOutlier = isOutsideOnTrackCut(id, residual, pull, 1);
            bool isRecoverable = m_recoverOutliers && !isOutlier && !isOutsideOnTrackCut(id, residual, pull, m_associationScaleFactor);
            bool flippedIsRecoverable =
                isMDT && flippedPull < pull - 0.1 && !isOutsideOnTrackCut(id, flippedResidual, flippedPull, m_associationScaleFactor);
            double innerRadius = 14.6;
            if (isMDT) {
                const MdtDriftCircleOnTrack* mdtdc = dynamic_cast<const MdtDriftCircleOnTrack*>(meas);
                if (mdtdc) innerRadius = mdtdc->detectorElement()->innerTubeRadius();
            }
 
            bool isDelta = isOutlier && isMDT && rTrackAbs < innerRadius && rTrackAbs > std::abs(rDrift);
 
            // remove all outliers that are too far from the track
            if (isOutlier) {
                if (isMDT) {
                    if (rTrackAbs > innerRadius) inBounds = false;
                } else if (pull > 10.) {
                    inBounds = false;
                }
            }
 
            std::unique_ptr<MdtDriftCircleOnTrack> mdtRotFlipped;
            std::unique_ptr<const CompetingMuonClustersOnTrack> updatedCompRot;
            bool flipSign = false;
            if (!pseudo) {
                const MdtDriftCircleOnTrack* mdtRot = isMDT ? dynamic_cast<const MdtDriftCircleOnTrack*>(meas) : nullptr;
                if (mdtRot && mdtRot->prepRawData() && mdtRot->prepRawData()->adc() < m_adcCut) {
                    isNoise = true;
                    isOutlier = true;
                    isRecoverable = false;
                    flippedIsRecoverable = false;
                    isDelta = false;
                }
                if (flippedIsRecoverable) {
                    residual = flippedResidual;
                    pull = flippedPull;
                    flipSign = true;
                }
                if (!isNoise && flipSign) {
                    if (mdtRot) {
                        mdtRotFlipped = std::make_unique<MdtDriftCircleOnTrack>(*mdtRot);
                        Trk::DriftCircleSide side = rDrift < 0. ? Trk::RIGHT : Trk::LEFT;
                        m_mdtRotCreator->updateSign(*mdtRotFlipped, side);
                        double rDriftFlip = mdtRotFlipped->localParameters()[Trk::locR];
                        int signRot = rDrift < 0 ? -1 : 1;
                        int signRotFlip = rDriftFlip < 0 ? -1 : 1;
                        if (rDrift != 0. && signRot == signRotFlip) {
                            ATH_MSG_WARNING(" failed to flip sign of MDT " << rDrift << "  flip " << rDriftFlip);
                        }
                    } else {
                        ATH_MSG_WARNING(" failed to dynamic_cast measurement with MDT Identifier to a MdtDriftCircleOnTrack");
                    }
                } else {
                    // if the outlier MDT is not recoverable by flipping the sign, add it as the worst outlier
                    if (isOutlier && isMDT && pull > largestmdtpull && !tsit->type(Trk::TrackStateOnSurface::Outlier)) {
                        largestmdtpull = pull;
                        mdtmeas = meas;
                    }
                }
 
                if (measuresPhi) {
                    bool isRpc = m_idHelperSvc->isRpc(id);
                    if (isRpc) {
                        int layer = 0;
                        StIndex stIndex = m_idHelperSvc->stationIndex(id);
                        if (stIndex == StIndex::BM && m_idHelperSvc->rpcIdHelper().doubletR(id) == 1)
                            layer = 1;
                        else if (stIndex == StIndex::BO)
                            layer = 2;
                        rpcLayers.insert(layer);
                    }
 
                    bool isTgc = m_idHelperSvc->isTgc(id);
                    if (isTgc) {
                        int layer = 0;
                        StIndex stIndex = m_idHelperSvc->stationIndex(id);
                        if (stIndex == StIndex::EM) {
                            std::string stName = m_idHelperSvc->chamberNameString(id);
                            if (stName[1] == '1')
                                layer = 1;
                            else if (stName[1] == '2')
                                layer = 2;
                            else if (stName[1] == '3')
                                layer = 3;
                            else {
                                ATH_MSG_WARNING("Unable to calculate TGC layer for " << m_idHelperSvc->toString(id));
                                layer = -1;
                            }
                        }
                        if (layer != -1) tgcLayers.insert(layer);
                    }
                    if (m_idHelperSvc->issTgc(id)) {
                        int layer = 4;
                        if (m_idHelperSvc->stgcIdHelper().multilayer(id) == 2) layer = 5;
                        tgcLayers.insert(layer);
                        ATH_MSG_VERBOSE("adding STGC phi hit " << layer << " size " << tgcLayers.size());
                    }
                }
 
                if (m_cleanCompROTs) {
                    const CompetingMuonClustersOnTrack* crot = (measuresPhi && !isMDT && m_idHelperSvc->isRpc(id))
                                                                   ? dynamic_cast<const CompetingMuonClustersOnTrack*>(meas)
                                                                   : nullptr;
                    if (crot) {
                        ATH_MSG_DEBUG(" CompetingMuonClustersOnTrack with rots " << crot->numberOfContainedROTs());
                        double minpos = 0.;
                        double minres = 0.;
                        double maxres = 0.;
                        double absminres = 0.;
                        double absmaxres = 0.;
                        for (unsigned int i = 0; i < crot->numberOfContainedROTs(); ++i) {
                            const MuonClusterOnTrack* cluster = &crot->rioOnTrack(i);
                            if (!cluster) continue;
                            double residual = cluster->localParameters()[Trk::locX] - pars->parameters()[Trk::locX];
                            double absres = residual < 0. ? -1. * residual : residual;
                            if (i == 0) {
                                minres = residual;
                                maxres = residual;
                                minpos = cluster->localParameters()[Trk::locX];
                                absminres = absres;
                                absmaxres = absres;
                            } else if (absres < absminres) {
                                minres = residual;
                                absminres = absres;
                                minpos = cluster->localParameters()[Trk::locX];
                            } else if (absres > absmaxres) {
                                maxres = residual;
                                absmaxres = absres;
                            }
                            ATH_MSG_VERBOSE(" ROT " << m_idHelperSvc->toString(cluster->identify()) << " lpos "
                                                    << cluster->localParameters()[Trk::locX] << " pars " << pars->parameters()[Trk::locX]
                                                    << " residual " << residual);
                        }
                        ATH_MSG_DEBUG(" residuals: min  " << minres << " max " << maxres << " diff " << maxres - minres);
                        bool splitCompRot = false;
                        if (std::abs(maxres - minres) > 100 && absmaxres - absminres > 20 && crot->numberOfContainedROTs() < 20) {
                            ATH_MSG_DEBUG(" recoverable double cluster ");
                            splitCompRot = true;
                        }
                        if (splitCompRot) {
                            std::list<const Trk::PrepRawData*> prdList;
                            ATH_MSG_DEBUG(" Splitting comp rot ");
                            for (unsigned int i = 0; i < crot->numberOfContainedROTs(); ++i) {
                                const MuonClusterOnTrack* cluster = &crot->rioOnTrack(i);
                                if (!cluster) continue;
                                double residual = cluster->localParameters()[Trk::locX] - minpos;
                                double absres = residual < 0. ? -1. * residual : residual;
                                if (absres < 40) {
                                    ATH_MSG_DEBUG("   NEW ROT " << m_idHelperSvc->toString(cluster->identify()) << " lpos "
                                                                << cluster->localParameters()[Trk::locX] << " pars "
                                                                << pars->parameters()[Trk::locX] << " residual " << residual);
                                    prdList.push_back(cluster->prepRawData());
                                }
                            }
                            if (prdList.empty()) {
                                ATH_MSG_WARNING("No clusters selected during comprot cleaning, keeping old cluster");
                            } else {
                                updatedCompRot=m_compRotCreator->createBroadCluster(prdList, 0.);
                                ++state.numberOfCleanedCompROTs;
                            }
                        }
                    }
                }
            }
 
            state.measInfo.emplace_back(id, chId, chIndex, inBounds, residual, pull, tsit, meas, pars, std::move(resPull), nullptr);
            MCTBCleaningInfo& info = state.measInfo.back();
            if (flipSign) { info.flippedMdt = std::move(mdtRotFlipped); }
            info.isNoise = isNoise;
            if (updatedCompRot) {
                ATH_MSG_DEBUG("updated competing ROT");
                info.cleanedCompROT = std::move(updatedCompRot);
                if (info.cleanedCompROT->associatedSurface() != meas->associatedSurface()) {
                     std::unique_ptr<Trk::TrackParameters> exPars = m_extrapolator->extrapolate(ctx,
                                                                        *pars, info.cleanedCompROT->associatedSurface(),
                                                                        Trk::anyDirection, false, Trk::muon);
                    if (!exPars) {
                        ATH_MSG_WARNING("Update of comp rot parameters failed, keeping old ones");
                        info.cleanedCompROT.reset();
                    } else {
                        info.pars = exPars.get();
                        state.parsToBeDeleted.emplace_back(std::move(exPars));
                    }
                }
            }
 
            ChamberLayerStatistics& chamberStatistics = state.chamberLayerStatistics[chIndex];
            chamberStatistics.chIndex = chIndex;
 
            if (tsit->type(Trk::TrackStateOnSurface::Outlier)) {
                ATH_MSG_DEBUG("  Outlier " << m_idHelperSvc->toString(id) << m_printer->print(*(info.resPull)));
                if (isMDT)
                    ATH_MSG_DEBUG(" r_drift " << rDrift << "  r_trk " << rTrack << " dr " << std::setw(7)
                                              << Amg::error(meas->localCovariance(), Trk::locR));
                if (isNoise) ATH_MSG_DEBUG(" noise ");
                if (!inBounds)
                    ATH_MSG_DEBUG(" outBounds ");
                else {
                    ATH_MSG_DEBUG(" inBounds  ");
                    if (flippedIsRecoverable)
                        ATH_MSG_DEBUG(" flipped is recoverable, rFlipped " << info.flippedMdt->localParameters()[Trk::locR]);
                    else if (isRecoverable)
                        ATH_MSG_DEBUG(" recoverable ");
                    else if (isDelta)
                        ATH_MSG_DEBUG(" delta ");
                }
 
                if (!inBounds) {
                    state.hasOfBoundsOutliers = true;
                    ++chamberStatistics.noutBounds;
                } else if (flippedIsRecoverable || isRecoverable)
                    ++chamberStatistics.nrecoverableOutliers;
                else if (isDelta)
                    ++chamberStatistics.ndeltas;
                else if (pull < 10)
                    ++chamberStatistics.noutliers;
 
                if (!pseudo) ++state.noutliers;  // don't count pseudos here
                info.useInFit = 0;
                continue;
            }
 
            // if we get here could be as hit
            ++chamberStatistics.nhits;
            if (!pseudo) ++state.nhits;  // pseudo's shouldn't count here
 
            if (flipSign) ++state.numberOfFlippedMdts;
 
            std::string idString = pseudo ? " Pseudomeasurement " : m_idHelperSvc->toString(id);
            std::string boundCheck = inBounds ? "inBounds" : "outBounds";
            ATH_MSG_VERBOSE(m_printer->print(*(info.resPull)) << " " << idString << " " << boundCheck);
            if (isNoise)
                ATH_MSG_VERBOSE(" noise");
            else if (isOutlier)
                ATH_MSG_VERBOSE(" outlier");
            if (isMDT) {
                ATH_MSG_VERBOSE(" r_drift " << std::setw(8) << rDrift << "  r_trk " << std::setw(7) << rTrack << " dr " << std::setw(7)
                                            << Amg::error(meas->localCovariance(), Trk::locR));
                if (flipSign)
                    ATH_MSG_VERBOSE(" flipped rDrift " << info.flippedMdt->localParameters()[Trk::locR] << "  pull " << flippedPull);
                else if (isDelta)
                    ATH_MSG_VERBOSE(" delta ");
            } else {
                const RpcClusterOnTrack* rpc = dynamic_cast<const RpcClusterOnTrack*>(meas);
                if (!rpc) {
                    const CompetingMuonClustersOnTrack* crot = dynamic_cast<const CompetingMuonClustersOnTrack*>(meas);
                    if (crot) { rpc = dynamic_cast<const RpcClusterOnTrack*>(crot->containedROTs().front()); }
                }
                if (rpc) ATH_MSG_DEBUG(" time " << rpc->prepRawData()->time() - rpc->globalPosition().mag() / 300.);
            }
 
            if (pseudo) {
                // pseudo measurements should not be included in chamber hit removal so stop here
                if (!isNoise && isOutlier) { state.largePullPseudoMeasurements.insert(meas); }
                continue;
 
            } else {
                // count hits
                EtaPhiHits& hits = state.hitsPerChamber[info.chId];
                if (measuresPhi)
                    ++hits.nphi;
                else
                    ++hits.neta;
 
                // bound checks
                if (!inBounds) {
                    EtaPhiHits& outHits = state.outBoundsPerChamber[info.chId];
                    if (measuresPhi)
                        ++outHits.nphi;
                    else
                        ++outHits.neta;
                }
 
                // measurements with large pull
                if (!tsit->type(Trk::TrackStateOnSurface::Outlier) && isOutlier && !isMDT) { state.largePullMeasurements.insert(meas); }
 
                // pulls per chamber
                if (!m_idHelperSvc->isTrigger(info.chId)) {
                    ChamberPullInfo& pullInfoHits = measuresPhi ? state.pullSumPhi : state.pullSum;
                    pullInfoHits.pullSum += pull;
                    ++pullInfoHits.nhits;
                    if (pull > pullInfoHits.maxPull) pullInfoHits.maxPull = pull;
 
                    ChamberPullInfo& pullInfoCh = measuresPhi ? state.pullSumPerChamberPhi[info.chId] : state.pullSumPerChamber[info.chId];
                    pullInfoCh.pullSum += pull;
                    ++pullInfoCh.nhits;
                    if (pull > pullInfoCh.maxPull) pullInfoCh.maxPull = pull;
                } else {
                    ChamberPullInfo& pullInfoTrig = measuresPhi ? state.pullSumPhi : state.pullSumTrigEta;
                    if (!measuresPhi) state.pullSumPerChamberEta[info.chId];
                    pullInfoTrig.pullSum += pull;
                    ++pullInfoTrig.nhits;
                    if (pull > pullInfoTrig.maxPull) pullInfoTrig.maxPull = pull;
                }
            }
        }
 
        // if there was an MDT outlier add it to the list of hits to be removed
        if (mdtmeas) state.largePullMeasurements.insert(mdtmeas);
 
        // check whether we have sufficient layers to savely clean the RPC/TGC comp rots
        unsigned int nphiLayers = rpcLayers.size() + tgcLayers.size();
        if (state.hasVertexConstraint) { nphiLayers += 1; }
        if (state.nIdHits > 0) { nphiLayers += 2; }
 
        if ((nphiLayers > 2) && state.numberOfCleanedCompROTs > 0) {
            ATH_MSG_DEBUG(" Sufficient constraints to clean competing ROTs: trigger layers " << nphiLayers);
        } else {
            // reset counter so no cleaning is performed
            state.numberOfCleanedCompROTs = 0;
        }
 
        // update sl fit configuration if track has ID hits or vertex constraint
        if (!m_use_slFit && state.slFit && (state.hasVertexConstraint || state.nIdHits > 0) && fieldCache.solenoidOn()) {
            state.slFit = false;
        }
 
        if (state.hasVertexConstraint || state.nIdHits) {
            if (state.hasVertexConstraint) ATH_MSG_DEBUG(" Track has vertex contraint:");
            if (state.nIdHits > 0) ATH_MSG_DEBUG(" Track has  ID Hits: " << state.nIdHits);
            if (fieldCache.solenoidOn())
                ATH_MSG_DEBUG(" Solenoid On");
            else
                ATH_MSG_DEBUG(" Solenoid Off");
            if (fieldCache.toroidOn())
                ATH_MSG_DEBUG(" Toroid On");
            else
                ATH_MSG_DEBUG(" Toroid Off");
            if (state.slFit)
                ATH_MSG_DEBUG(" Use SL fit");
            else
                ATH_MSG_DEBUG(" Use curved fit");
        }
 
        std::set<Identifier> chambersToBeRemoved, chambersToBeRemovedPhi, goodEtaLayers, goodPhiLayers;
 
        // check whether there are any chambers in eta that should be removed
        extractChambersToBeRemoved(state, chambersToBeRemoved);
 
        // check whether there are any chambers in phi that should be removed
        extractChambersToBeRemoved(state, chambersToBeRemovedPhi, true);
 
        ATH_MSG_DEBUG("Chambers on track summary:");
 
        EtaPhiPerChamberIt chit = state.hitsPerChamber.begin();
        EtaPhiPerChamberIt chit_end = state.hitsPerChamber.end();
        for (; chit != chit_end; ++chit) {
            if (!chit->first.is_valid()) continue;
 
            bool isPrec = m_idHelperSvc->isMdt(chit->first) || m_idHelperSvc->isCsc(chit->first) || m_idHelperSvc->isMM(chit->first) ||
                          m_idHelperSvc->issTgc(chit->first);
            if (isPrec) {
                StIndex stIndex = m_idHelperSvc->stationIndex(chit->first);
                state.stations.insert(stIndex);
            }
 
            EtaPhiHits& nhits = chit->second;
            EtaPhiHits& noutBounds = state.outBoundsPerChamber[chit->first];
            ATH_MSG_DEBUG(" chamber " << m_idHelperSvc->toStringChamber(chit->first) << "  eta hits " << nhits.neta << " outBounds "
                                      << noutBounds.neta << "  phi hits " << nhits.nphi << " outBounds " << noutBounds.nphi);
            if (nhits.neta != 0 && nhits.neta == noutBounds.neta) {
                ATH_MSG_DEBUG("   --> All eta hits out of bounds ");
                // remove eta hits for this chamber
                double fakePull = 1000. * nhits.neta;
                if (!chambersToBeRemoved.count(chit->first)) {
                    // just add it
                    state.chambersToBeRemoved.emplace_back(fakePull, chit->first);
                    chambersToBeRemoved.insert(chit->first);
                } else {
                    // replace existing entry if new one has larger Pull
                    PullChIt pIt = state.chambersToBeRemoved.begin();
                    PullChIt pIt_end = state.chambersToBeRemoved.end();
                    for (; pIt != pIt_end; ++pIt) {
                        if (pIt->second == chit->first && pIt->first < fakePull) {
                            pIt->first = fakePull;
                            break;
                        }
                    }
                }
            } else {
                if (noutBounds.neta != 0) ATH_MSG_DEBUG("   --> Some eta hits out of bounds ");
 
                if (isPrec && nhits.neta > 0) {
                    if (m_idHelperSvc->isSmallChamber(chit->first))
                        state.hasSmall = true;
                    else
                        state.hasLarge = true;
                }
            }
            if (nhits.nphi != 0 && nhits.nphi == noutBounds.nphi) {
                ATH_MSG_DEBUG("   --> All phi hits out of bounds ");
                // remove phi hits for this chamber
                double fakePull = 1000. * nhits.nphi;
                if (!chambersToBeRemovedPhi.count(chit->first)) {
                    // just add it
                    state.chambersToBeRemovedPhi.emplace_back(fakePull, chit->first);
                    chambersToBeRemovedPhi.insert(chit->first);
                } else {
                    // replace existing entry if new one has larger pull
                    PullChIt pIt = state.chambersToBeRemovedPhi.begin();
                    PullChIt pIt_end = state.chambersToBeRemovedPhi.end();
                    for (; pIt != pIt_end; ++pIt) {
                        if (pIt->second == chit->first) {
                            pIt->first = std::max(pIt->first, fakePull);
                            break;
                        }
                    }
                }
            } else {
                if (noutBounds.nphi != 0) ATH_MSG_DEBUG("   --> Some phi hits out of bounds ");
                if (nhits.nphi > 0) {
                    state.phiLayers.insert(m_idHelperSvc->phiIndex(chit->first));
                }
            }
        }
 
        state.nPhiConstraints = state.phiLayers.size() + state.nPseudoMeasurements;
        if (state.hasLarge && state.hasSmall) ++state.nPhiConstraints;
        if (state.stations.size() == 1 && nphiLayers > 1) ++state.nPhiConstraints;
 
        double pull_precisionhits = -1.;
        if (state.pullSum.nhits != 0) { pull_precisionhits = state.pullSum.pullSum / state.pullSum.nhits; }
        double pull_triggerhits = -1.;
        if (state.pullSumTrigEta.nhits != 0) { pull_triggerhits = state.pullSumTrigEta.pullSum / state.pullSumTrigEta.nhits; }
        double pull_phihits = -1.;
        if (state.pullSumPhi.nhits != 0) { pull_phihits = state.pullSumPhi.pullSum / state.pullSumPhi.nhits; }
        ATH_MSG_DEBUG("    pull sum: precision hits "
                      << pull_precisionhits << " trigger eta " << pull_triggerhits << "  phi " << pull_phihits << "  measurements "
                      << state.measInfo.size() << std::endl
                      << "  precision chambers " << state.pullSumPerChamber.size() << "  hits with large pull "
                      << state.largePullMeasurements.size() << "  pseudos with large pull " << state.largePullPseudoMeasurements.size()
                      << "  chambers to be removed " << state.chambersToBeRemoved.size() << "  phi " << state.chambersToBeRemovedPhi.size()
                      << "  phi lay " << state.phiLayers.size() << "  phi constr " << state.nPhiConstraints);
 
        if (state.pullSumPhi.nhits != 0) {
            double phiPull = state.pullSumPhi.pullSum / state.pullSumPhi.nhits;
            if (state.stations.size() == 1 && phiPull > 5.) { ATH_MSG_DEBUG(" Single station track with large phi pulls!! "); }
        } else
            ATH_MSG_DEBUG(" Track with no phi hits!! ");
    }

initialize()

StatusCode Muon::MuonTrackCleaner::initialize ( )

override

AlgTool initialize.

Definition at line 37 of file MuonTrackCleaner.cxx.

                                            {
        ATH_CHECK(m_trackFitter.retrieve());
        ATH_CHECK(m_slTrackFitter.retrieve());
        ATH_CHECK(m_idHelperSvc.retrieve());
        ATH_CHECK(m_edmHelperSvc.retrieve());
        ATH_CHECK(m_printer.retrieve());
        ATH_CHECK(m_extrapolator.retrieve());
        ATH_CHECK(m_pullCalculator.retrieve());
        ATH_CHECK(m_mdtRotCreator.retrieve());
        ATH_CHECK(m_compRotCreator.retrieve());
        ATH_CHECK(m_measurementUpdator.retrieve());
        ATH_CHECK(m_fieldCacheCondObjInputKey.initialize());
 
        return StatusCode::SUCCESS;
    }

isOutsideOnTrackCut()

bool Muon::MuonTrackCleaner::isOutsideOnTrackCut ( const Identifier &  id, double  res, double  pull, double  cutScaleFactor  ) const

private

check whether hit is an outlier

Definition at line 1473 of file MuonTrackCleaner.cxx.

                                                                                                                         {
        bool isMdt = id.is_valid() ? m_idHelperSvc->isMdt(id) : false;
        bool measuresPhi = id.is_valid() ? m_idHelperSvc->measuresPhi(id) : false;
 
        // look at pulls and identify outliers
        double pullCut = measuresPhi ? m_pullCutPhi : m_pullCut;
 
        if (isMdt) {
            // if mdt residual cut is activated check whether the residual is small that 80% of the cut of
            if (m_useMdtResiCut) {
                return std::abs(res) > cutScaleFactor * m_mdtResiCut;
            } else {
                return std::abs(pull) > cutScaleFactor * pullCut;
            }
        } else {
            return std::abs(pull) > cutScaleFactor * pullCut;
        }
    }

outlierRecovery()

std::unique_ptr< Trk::Track > Muon::MuonTrackCleaner::outlierRecovery ( const EventContext &  ctx, std::unique_ptr< Trk::Track >  track, CleaningState &  state, const MuonStationIndex::ChIndex *  currentIndex = nullptr  ) const

private

recover outliers that are within the cuts.

If the additional chamber index is provided, the code will only consider that particular layer

Definition at line 656 of file MuonTrackCleaner.cxx.

                                                                                                                     {
        const Trk::Perigee* perigee = track->perigeeParameters();
        if (!perigee) {
            ATH_MSG_DEBUG("   track without perigee ");
            return nullptr;
        }
 
        ATH_MSG_VERBOSE("  outlierRecovery: ");
        if (currentIndex) ATH_MSG_VERBOSE(" layer " << chName(*currentIndex));
        ATH_MSG_VERBOSE("  printing chamber statistics ");
 
        std::set<ChIndex> recoverableLayers;
        for (auto [stationIndex, layer] : state.chamberLayerStatistics) {
            if (stationIndex == ChIndex::ChUnknown) continue;
            ATH_MSG_VERBOSE(print(layer));
 
            // skip all chamber layers except the requested one
            if (currentIndex && *currentIndex != stationIndex) continue;
 
            ChamberLayerStatistics& statistics = layer;
            unsigned int nhits = statistics.nhits;
            unsigned int noutliers = statistics.noutliers;
            // unsigned int ndeltas = statistics.ndeltas;
            unsigned int nrecoverableOutliers = statistics.nrecoverableOutliers;
            unsigned int noutBounds = statistics.noutBounds;
 
            if (nrecoverableOutliers > 0) {
                if (nhits + nrecoverableOutliers > 2 &&
                    ((noutBounds == 0 && noutliers == 0) || (nrecoverableOutliers != 0 && noutliers < 2))) {
                    recoverableLayers.insert(stationIndex);
                    ATH_MSG_DEBUG("   found recoverable layer " << chName(statistics.chIndex));
                }
            }
        }
 
        // if there is nothing to be done exit
        if (recoverableLayers.empty() && !state.hasOfBoundsOutliers) return track;
        bool addedHits = false;
        unsigned int removedOutOfBoundsHits(0);
 
        auto tsos = std::make_unique<Trk::TrackStates>();
        tsos->reserve(state.measInfo.size());
 
        // loop over hits
        for (auto& hit : state.measInfo) {
            if (!hit.useInFit) {
                if (hit.inBounds) {
                    if (recoverableLayers.count(hit.chIndex)) {
                        // check whether we can savely add hits in this chamber to the track
                        bool recover = !isOutsideOnTrackCut(hit.id, hit.residual, hit.pull, m_associationScaleFactor);
                        if (recover && m_onlyUseHitErrorInRecovery && hit.pars) {
                            std::optional<const Trk::ResidualPull> resPull{
                                m_pullCalculator->residualPull(hit.meas, hit.pars, Trk::ResidualPull::HitOnly)};
                            if (!resPull) {
                                ATH_MSG_DEBUG(" calculation of residual/pull failed !!!!! ");
                                recover = false;
                            } else {
                                recover = !isOutsideOnTrackCut(hit.id, resPull->residual().front(), std::abs(resPull->pull().front()),
 
                                                               m_associationScaleFactor);
                            }
                        }
                        if (recover) {
                            ATH_MSG_DEBUG("   adding outlier " << m_idHelperSvc->toString(hit.id) << " pull " << std::setw(7) << hit.pull);
                            if (hit.flippedMdt) {
                                double rDrift = hit.meas->localParameters()[Trk::locR];
                                double rDriftFlip = hit.flippedMdt->localParameters()[Trk::locR];
                                double rTrack = hit.pars->parameters()[Trk::locR];
                                ATH_MSG_DEBUG(" flipped MDT: r_orig " << rDrift << " flip " << rDriftFlip << " rTrack " << rTrack);
                            }
                            addedHits = true;
                            const Trk::MeasurementBase* newMeas = hit.flippedMdt ? hit.flippedMdt.get() : hit.meas;
 
                            tsos->push_back(MuonTSOSHelper::cloneTSOSWithUpdate(*hit.originalState, *newMeas, *hit.pars,
                                                                               Trk::TrackStateOnSurface::Measurement));
                        } else {
                            tsos->push_back(hit.originalState->clone());
                        }
                    }
                    // layer not recoverable, drop the outliers: but if RPC, TGC, or CSC, expect track to go through all layers, so add a
                    // hole instead
                    if (m_idHelperSvc->isRpc(hit.id) || m_idHelperSvc->isTgc(hit.id) || m_idHelperSvc->isCsc(hit.id))
                        tsos->push_back(MuonTSOSHelper::createHoleTSOS(hit.pars->uniqueClone()));
                } else {
                    ++removedOutOfBoundsHits;
                    // if RPC, TGC, or CSC, expect track to go through all layers: add a hole to replace lost outlier
                    if (m_idHelperSvc->isRpc(hit.id) || m_idHelperSvc->isTgc(hit.id) || m_idHelperSvc->isCsc(hit.id))
                        tsos->push_back(MuonTSOSHelper::createHoleTSOS(hit.pars->uniqueClone()));
                    ATH_MSG_DEBUG("   removing out of bounds outlier " << m_idHelperSvc->toString(hit.id) << " pull " << std::setw(7)
                                                                       << hit.pull);
                }
            } else {
                tsos->push_back(hit.originalState->clone());
            }
        }
 
        if (!addedHits && removedOutOfBoundsHits == 0) {
            ATH_MSG_DEBUG(" track unchanged ");
            return track;
        }
 
        if (tsos->size() < 6) {
            ATH_MSG_WARNING(" too few hits, cannot add hits. This should not happen ");
            return nullptr;
        }
 
        // create new track
        std::unique_ptr<Trk::Track> cleanedTrack =
            std::make_unique<Trk::Track>(track->info(), std::move(tsos), track->fitQuality() ? track->fitQuality()->uniqueClone() : nullptr);
 
        if (!addedHits) {
            ATH_MSG_DEBUG(" only removed out of bound hits, returning track without new fit ");
            return cleanedTrack;
        }
 
        // fit new track
        printStates(cleanedTrack.get());
 
        std::unique_ptr<Trk::Track> newTrack = fitTrack(ctx, *cleanedTrack, track->info().particleHypothesis(), state.slFit);
 
        if (newTrack) {
            ATH_MSG_DEBUG("Outlier recovery successfull ");
            init(ctx, *newTrack, state);
        } else {
            ATH_MSG_DEBUG("refit after outlier recovery failed ");
        }
 
        return newTrack;
    }

print()

std::string Muon::MuonTrackCleaner::print ( MuonTrackCleaner::ChamberLayerStatistics &  statistics )

staticprivate

Definition at line 1492 of file MuonTrackCleaner.cxx.

                                                                                        {
        std::ostringstream sout;
 
        unsigned int nhits = statistics.nhits;
        unsigned int noutliers = statistics.noutliers;
        unsigned int ndeltas = statistics.ndeltas;
        unsigned int nrecoverableOutliers = statistics.nrecoverableOutliers;
        unsigned int noutBounds = statistics.noutBounds;
 
        sout << chName(statistics.chIndex) << " hits " << std::setw(6) << nhits << " outliers " << std::setw(6)
             << noutliers << " deltas " << std::setw(6) << ndeltas << " recoverableOutliers " << std::setw(6) << nrecoverableOutliers
             << " outBounds " << std::setw(6) << noutBounds;
 
        return sout.str();
    }

printStates()

void Muon::MuonTrackCleaner::printStates ( Trk::Track *  track ) const

private

Definition at line 1542 of file MuonTrackCleaner.cxx.

                                                            {
        const Trk::TrackStates* states = track->trackStateOnSurfaces();
        int nStates = 0;
        if (states) nStates = states->size();
        ATH_MSG_DEBUG("Calling fit with hits: " << nStates<<std::endl<<m_printer->printMeasurements(*track));
    }

recoverFlippedMdt()

std::unique_ptr< Trk::Track > Muon::MuonTrackCleaner::recoverFlippedMdt ( const EventContext &  ctx, std::unique_ptr< Trk::Track >  track, CleaningState &  state  ) const

private

flip signs of MDT hits with large pull if pull if the oppositely signed radius is small

Definition at line 235 of file MuonTrackCleaner.cxx.

                                                                                                {
        if (!m_flipMdtDriftRadii || state.numberOfFlippedMdts == 0) return track;
 
        const Trk::Perigee* perigee = track->perigeeParameters();
        if (!perigee) {
            ATH_MSG_DEBUG("   track without perigee ");
            return nullptr;
        }
 
        ATH_MSG_DEBUG(" Trying to flip MDT signs: total number of hits with wrong sign " << state.numberOfFlippedMdts);
 
        auto tsos = std::make_unique<Trk::TrackStates>();
        tsos->reserve(state.measInfo.size());
 
        unsigned int nmeas = 0;
        // loop over hits
        InfoIt hit = state.measInfo.begin();
        InfoIt hit_end = state.measInfo.end();
        for (; hit != hit_end; ++hit) {
            // hits that are flagged as outlier or hits in the chamber to be removed are added as Outlier
            if (!hit->useInFit) {
                ATH_MSG_DEBUG("   removing hit " << m_idHelperSvc->toString(hit->id) << " pull " << hit->resPull->pull().front());
                if (hit->inBounds)
                    tsos->push_back(MuonTSOSHelper::cloneTSOSWithUpdate(*hit->originalState, *hit->meas, *hit->pars,
                                                                       Trk::TrackStateOnSurface::Outlier));
 
                continue;
            } else {
                if (hit->meas) ++nmeas;
 
                if (hit->flippedMdt) {
                    tsos->push_back(MuonTSOSHelper::cloneTSOSWithUpdate(*hit->originalState, *hit->flippedMdt, *hit->pars,
                                                                       Trk::TrackStateOnSurface::Measurement));
                    ATH_MSG_DEBUG("   flipping sign hit " << m_idHelperSvc->toString(hit->id) << " pull " << hit->resPull->pull().front());
                } else {
                    tsos->push_back(hit->originalState->clone());
                }
            }
        }
 
        if (nmeas < 6) {
            ATH_MSG_DEBUG(" too few hits, cannot flip MDT hit sign ");
            return nullptr;
        }
 
        // create new track
        std::unique_ptr<Trk::Track> cleanedTrack =
            std::make_unique<Trk::Track>(track->info(), std::move(tsos), track->fitQuality() ? track->fitQuality()->uniqueClone() : nullptr);
        printStates(cleanedTrack.get());
 
        // fit new track
        std::unique_ptr<Trk::Track> newTrack = fitTrack(ctx, *cleanedTrack, track->info().particleHypothesis(), state.slFit);
 
        if (newTrack) {
            init(ctx, *newTrack, state);
            return newTrack;
        } else
            return nullptr;
    }

removeChamber()

MuonTrackCleaner::ChamberRemovalOutput Muon::MuonTrackCleaner::removeChamber ( const EventContext &  ctx, const std::unique_ptr< Trk::Track > &  track, Identifier  chId, bool  removePhi, bool  removeEta, CleaningState &  state  ) const

private

remove chamber from track

Definition at line 579 of file MuonTrackCleaner.cxx.

                                                                                                                                       {
        ATH_MSG_DEBUG(" removing chamber " << m_idHelperSvc->toStringChamber(chId));
 
        // store result
        ChamberRemovalOutput result;
        result.chId = chId;
        const Trk::Perigee* perigee = track->perigeeParameters();
        if (!perigee) {
            ATH_MSG_DEBUG("   track without perigee ");
            return result;
        }
 
        auto tsos = std::make_unique<Trk::TrackStates>();
        tsos->reserve(state.measInfo.size());
 
        unsigned int nmeas = 0;
        // loop over hits
        std::set<Identifier> stations{};
        for (MCTBCleaningInfo& hit : state.measInfo) {
            if (hit.id.is_valid()) {
                bool measuresPhi = m_idHelperSvc->measuresPhi(hit.id);
                bool remove = hit.chId == chId && ((removePhi && measuresPhi) || (removeEta && !measuresPhi));
                // hits that are flagged as outlier or hits in the chamber to be removed are added as Outlier
                if (!hit.useInFit || remove) {
                    ATH_MSG_DEBUG("   removing hit " << m_idHelperSvc->toString(hit.id) << " pull " << hit.resPull->pull().front());
                    // add as outlier
                    if (hit.inBounds)
                        tsos->push_back(MuonTSOSHelper::cloneTSOSWithUpdate(*hit.originalState, *hit.meas, *hit.pars,
                                                                           Trk::TrackStateOnSurface::Outlier));
 
                    // if removed, add hit to vector of hits
                    // but only if the hit was not already an outlier to be skipped!
                    if (remove && hit.useInFit) result.removedHits.push_back(&hit);
                    hit.useInFit = 0;
                    continue;
                }
            }
            if (hit.meas) {
                ++nmeas;
                // Through how many detector layers went actually the track ?
                if (m_idHelperSvc->isMdt(hit.id)) stations.insert(m_idHelperSvc->mdtIdHelper().multilayerID(hit.id));
                // Keep it for the run 2 legacy
                else if (m_idHelperSvc->isCsc(hit.id)) stations.insert(m_idHelperSvc->cscIdHelper().elementID(hit.id));
                // Multilayer defines here whether the track went through the IP or HO wedge
                else if (m_idHelperSvc->isMM(hit.id)) stations.insert(m_idHelperSvc->mmIdHelper().multilayerID(hit.id));
                // Hits on track in the site before or after the Micromega wedges... yeah micromega sandwich
                else if (m_idHelperSvc->issTgc(hit.id)) stations.insert(m_idHelperSvc->stgcIdHelper().multilayerID(hit.id));
            }
            tsos->push_back(hit.originalState->clone());
        }
 
        if (nmeas < 6 || stations.size() < 2) {
            ATH_MSG_DEBUG(" too few hits, cannot perform chamberCleaning ");
            return result;
        }
 
        // create new track
        std::unique_ptr<Trk::Track> cleanedTrack =
            std::make_unique<Trk::Track>(track->info(), std::move(tsos), track->fitQuality() ? track->fitQuality()->uniqueClone() : nullptr);
 
        // fit new track
        printStates(cleanedTrack.get());
 
        if (!cleanedTrack->perigeeParameters()) { ATH_MSG_DEBUG("   track without perigee "); }
 
        std::unique_ptr<Trk::Track> newTrack = fitTrack(ctx, *cleanedTrack, track->info().particleHypothesis(), state.slFit);
 
        if (newTrack && !newTrack->perigeeParameters()) {
            ATH_MSG_DEBUG("   fitted track without perigee " << *newTrack << " " << newTrack->perigeeParameters());
        }
 
        result.track = std::move(newTrack);
        return result;
    }

unremoveHits()

void Muon::MuonTrackCleaner::unremoveHits ( ChamberRemovalOutput &  result )

staticprivate

Definition at line 1537 of file MuonTrackCleaner.cxx.

                                                                    {
        // loop over removed hits and 'unremove' them so they are used in the next iteration
        for (auto* hit : result.removedHits) hit->useInFit = 1;
    }

Member Data Documentation

m_adcCut

Gaudi::Property<double> Muon::MuonTrackCleaner::m_adcCut {this, "AdcCut", 50.}

private

Definition at line 313 of file MuonTrackCleaner.h.

m_associationScaleFactor

Gaudi::Property<double> Muon::MuonTrackCleaner::m_associationScaleFactor {this, "AssociationScaleFactor", 0.7}

private

Definition at line 306 of file MuonTrackCleaner.h.

m_avePullSumPerChamberCut

Gaudi::Property<double> Muon::MuonTrackCleaner::m_avePullSumPerChamberCut {this, "MaxAvePullSumPerChamber", 3.5}

private

Definition at line 308 of file MuonTrackCleaner.h.

m_chi2Cut

Gaudi::Property<double> Muon::MuonTrackCleaner::m_chi2Cut {this, "Chi2Cut", 100.}

private

Definition at line 302 of file MuonTrackCleaner.h.

m_cleanCompROTs

Gaudi::Property<bool> Muon::MuonTrackCleaner::m_cleanCompROTs {this, "CleanCompROTs", true}

private

Definition at line 311 of file MuonTrackCleaner.h.

m_compRotCreator

ToolHandle<IMuonCompetingClustersOnTrackCreator> Muon::MuonTrackCleaner::m_compRotCreator

private

Initial value:

{
            this, "CompRotCreator", "Muon::TriggerChamberClusterOnTrackCreator/TriggerChamberClusterOnTrackCreator"}

Definition at line 288 of file MuonTrackCleaner.h.

m_edmHelperSvc

ServiceHandle<Muon::IMuonEDMHelperSvc> Muon::MuonTrackCleaner::m_edmHelperSvc

private

Initial value:

{this, "edmHelper", "Muon::MuonEDMHelperSvc/MuonEDMHelperSvc",
                                                              "Handle to the service providing the IMuonEDMHelperSvc interface"}

Definition at line 292 of file MuonTrackCleaner.h.

m_extrapolator

ToolHandle<Trk::IExtrapolator> Muon::MuonTrackCleaner::m_extrapolator {this, "Extrapolator", "Trk::Extrapolator/AtlasExtrapolator"}

private

Definition at line 296 of file MuonTrackCleaner.h.

m_fieldCacheCondObjInputKey

SG::ReadCondHandleKey<AtlasFieldCacheCondObj> Muon::MuonTrackCleaner::m_fieldCacheCondObjInputKey

private

Initial value:

{this, "AtlasFieldCacheCondObj", "fieldCondObj",
                                                                                  "Name of the Magnetic Field conditions object key"}

Definition at line 298 of file MuonTrackCleaner.h.

m_flipMdtDriftRadii

Gaudi::Property<bool> Muon::MuonTrackCleaner::m_flipMdtDriftRadii {this, "FlipMdtDriftRadii", true}

private

Definition at line 310 of file MuonTrackCleaner.h.

m_idHelperSvc

ServiceHandle<Muon::IMuonIdHelperSvc> Muon::MuonTrackCleaner::m_idHelperSvc {this, "MuonIdHelperSvc", "Muon::MuonIdHelperSvc/MuonIdHelperSvc"}

private

Definition at line 295 of file MuonTrackCleaner.h.

m_iterate

Gaudi::Property<bool> Muon::MuonTrackCleaner::m_iterate {this, "Iterate", 0.7}

private

Definition at line 314 of file MuonTrackCleaner.h.

m_mdtResiCut

Gaudi::Property<double> Muon::MuonTrackCleaner::m_mdtResiCut {this, "MdtResiCut", 1.}

private

Definition at line 304 of file MuonTrackCleaner.h.

m_mdtRotCreator

ToolHandle<Muon::IMdtDriftCircleOnTrackCreator> Muon::MuonTrackCleaner::m_mdtRotCreator

private

Initial value:

{this, "MdtRotCreator",
                                                                        "Muon::MdtDriftCircleOnTrackCreator/MdtDriftCircleOnTrackCreator"}

Definition at line 286 of file MuonTrackCleaner.h.

m_measurementUpdator

ToolHandle<Trk::IUpdator> Muon::MuonTrackCleaner::m_measurementUpdator {this, "MeasurementUpdator", "Trk::KalmanUpdator/MuonMeasUpdator"}

private

Definition at line 285 of file MuonTrackCleaner.h.

m_ncycles

Gaudi::Property<unsigned int> Muon::MuonTrackCleaner::m_ncycles {this, "CleaningCycles", 5}

private

Definition at line 307 of file MuonTrackCleaner.h.

m_onlyUseHitErrorInRecovery

Gaudi::Property<bool> Muon::MuonTrackCleaner::m_onlyUseHitErrorInRecovery {this, "OnlyUseHitErrorInRecovery", true}

private

Definition at line 312 of file MuonTrackCleaner.h.

m_printer

PublicToolHandle<Muon::MuonEDMPrinterTool> Muon::MuonTrackCleaner::m_printer {this, "Printer", "Muon::MuonEDMPrinterTool/MuonEDMPrinterTool"}

private

Definition at line 294 of file MuonTrackCleaner.h.

m_pullCalculator

ToolHandle<Trk::IResidualPullCalculator> Muon::MuonTrackCleaner::m_pullCalculator

private

Initial value:

{this, "PullCalculator",
                                                                  "Trk::ResidualPullCalculator/ResidualPullCalculator"}

Definition at line 290 of file MuonTrackCleaner.h.

m_pullCut

Gaudi::Property<double> Muon::MuonTrackCleaner::m_pullCut {this, "PullCut", 5.}

private

Definition at line 303 of file MuonTrackCleaner.h.

m_pullCutPhi

Gaudi::Property<double> Muon::MuonTrackCleaner::m_pullCutPhi {this, "PullCutPhi", 10.}

private

Definition at line 305 of file MuonTrackCleaner.h.

m_recoverOutliers

Gaudi::Property<bool> Muon::MuonTrackCleaner::m_recoverOutliers {this, "RecoverOutliers", true}

private

Definition at line 309 of file MuonTrackCleaner.h.

m_slTrackFitter

ToolHandle<Trk::ITrackFitter> Muon::MuonTrackCleaner::m_slTrackFitter {this, "SLFitter", "Trk::GlobalChi2Fitter/MCTBSLFitterMaterialFromTrack"}

private

Definition at line 284 of file MuonTrackCleaner.h.

m_trackFitter

ToolHandle<Trk::ITrackFitter> Muon::MuonTrackCleaner::m_trackFitter

private

Initial value:

{
            this,
            "Fitter",
            "Trk::GlobalChi2Fitter/MCTBFitterMaterialFromTrack",
        }

Definition at line 279 of file MuonTrackCleaner.h.

m_use_slFit

Gaudi::Property<bool> Muon::MuonTrackCleaner::m_use_slFit {this, "UseSLFit", false}

private

Definition at line 315 of file MuonTrackCleaner.h.

m_useMdtResiCut

Gaudi::Property<bool> Muon::MuonTrackCleaner::m_useMdtResiCut {this, "UseMdtResiCut", false}

private

Definition at line 301 of file MuonTrackCleaner.h.

---

The documentation for this class was generated from the following files:

• MuonTrackCleaner.h
• MuonTrackCleaner.cxx