The SegmentLineFitter is a standalone module to fit a straight line to calibrated muon space points. More...

#include <SegmentLineFitter.h>

Inheritance diagram for MuonR4::SegmentFit::SegmentLineFitter:

Collaboration diagram for MuonR4::SegmentFit::SegmentLineFitter:

Classes
struct	ConfigSwitches
	Configuration object of the ATLAS implementation. More...
struct	Config
	Full configuration object. More...

Public Types
using	Fitter_t = Acts::Experimental::CompositeSpacePointLineFitter
	Abrivation of the actual line fitter.
using	LinePar_t = Fitter_t::ParamVec_t
	Abrivation of the fitted line parameters.
using	Hit_t = Segment::MeasType
	Abrivation of the space point type to use.
using	HitVec_t = std::vector<Hit_t>
	Collection of space points.
using	FitPars_t = Fitter_t::FitParameters
	Abrivation of the fit parameters.
using	Result_t = Fitter_t::FitResult<HitVec_t>
	Abrivation of the fit result.
using	FitOpts_t = Fitter_t::FitOptions<HitVec_t, ISpacePointCalibrator>
	Abrivation of the fit options.
using	Selector_t = Fitter_t::Selector_t<CalibratedSpacePoint>
	Abrivation of the hit selector to choose valid hits.
using	HitState = CalibratedSpacePoint::State
	Abrivation of the fit state flag.

Public Member Functions
	SegmentLineFitter (const std::string &name, Config &&config)
	Standard constructor.
std::unique_ptr< Segment >	fitSegment (const EventContext &ctx, const SegmentSeed *parent, const LinePar_t &startPars, const Amg::Transform3D &localToGlobal, HitVec_t &&calibHits) const
	Fit a set of measurements to a straight segment line.
bool	msgLvl (const MSG::Level lvl) const
	Test the output level.
MsgStream &	msg () const
	The standard message stream.
MsgStream &	msg (const MSG::Level lvl) const
	The standard message stream.
void	setLevel (MSG::Level lvl)
	Change the current logging level.

Private Member Functions
Result_t	callLineFit (const Acts::CalibrationContext &cctx, const Parameters &startPars, const Amg::Transform3D &localToGlobal, HitVec_t &&calibHits) const
	Calls the underlying line fitter to determine the segment parameters.
bool	removeOutliers (const Acts::CalibrationContext &cctx, const SegmentSeed &seed, const Amg::Transform3D &localToGlobal, const LinePar_t &startPars, Result_t &fitResult) const
	Cleans the fitted segment from the most outlier hit and then attempts to refit the segment.
bool	plugHoles (const Acts::CalibrationContext &cctx, const SegmentSeed &seed, const Amg::Transform3D &localToGlobal, Result_t &toRecover) const
	Recovery of missed hits.
void	eraseWrongHits (Result_t &candidate) const
	Removes all hits from the segment which are obvious outliers.
void	cleanStripLayers (HitVec_t &hits) const
	Marks duplicate hits on a strip layer as outliers to avoid competing contributions from the same layers in the fit.
bool	betterResult (const Result_t &newResult, const Result_t &oldResult) const
	Returns whether the new fit result is better than the one from the previous iteration.
std::unique_ptr< Segment >	convertToSegment (const Amg::Transform3D &locToGlobTrf, const SegmentSeed *parentSeed, Result_t &&toConvert) const
	Converts the fit result into a segment object.
bool	checkPrecHitCount (const HitVec_t &candidateHits) const
	Checks if the candidate has enough precision hits to fit a segment.
void	initMessaging () const
	Initialize our message level and MessageSvc.

Private Attributes
Fitter_t	m_fitter
	Actual implementation of the straight line fit.
ConfigSwitches	m_cfg {}
	Configuration switches of the ATLAS fitter implementation.
Selector_t	m_goodHitSel {}
	Selector to identify the valid hits.
std::string	m_nm
	Message source name.
boost::thread_specific_ptr< MsgStream >	m_msg_tls
	MsgStream instance (a std::cout like with print-out levels).
std::atomic< IMessageSvc * >	m_imsg { nullptr }
	MessageSvc pointer.
std::atomic< MSG::Level >	m_lvl { MSG::NIL }
	Current logging level.
std::atomic_flag m_initialized	ATLAS_THREAD_SAFE = ATOMIC_FLAG_INIT
	Messaging initialized (initMessaging).

Detailed Description

The SegmentLineFitter is a standalone module to fit a straight line to calibrated muon space points.

The CompositeSpacePointLineFitter from the ACTS toolkit is used to perform the actual fit to the measurements. The SegmentLineFitter is a wrapper class taking care of relaunching fits of poor quality but with cleaned measurements and also to put back meaurements on the line that have been missed by the initial line fit.

Definition at line 28 of file SegmentLineFitter.h.

Member Typedef Documentation

◆ FitOpts_t

using MuonR4::SegmentFit::SegmentLineFitter::FitOpts_t = Fitter_t::FitOptions<HitVec_t, ISpacePointCalibrator>

Abrivation of the fit options.

Definition at line 43 of file SegmentLineFitter.h.

◆ FitPars_t

using MuonR4::SegmentFit::SegmentLineFitter::FitPars_t = Fitter_t::FitParameters

Abrivation of the fit parameters.

Definition at line 39 of file SegmentLineFitter.h.

◆ Fitter_t

using MuonR4::SegmentFit::SegmentLineFitter::Fitter_t = Acts::Experimental::CompositeSpacePointLineFitter

Abrivation of the actual line fitter.

Definition at line 31 of file SegmentLineFitter.h.

◆ Hit_t

using MuonR4::SegmentFit::SegmentLineFitter::Hit_t = Segment::MeasType

Abrivation of the space point type to use.

Definition at line 35 of file SegmentLineFitter.h.

◆ HitState

using MuonR4::SegmentFit::SegmentLineFitter::HitState = CalibratedSpacePoint::State

Abrivation of the fit state flag.

Definition at line 47 of file SegmentLineFitter.h.

◆ HitVec_t

using MuonR4::SegmentFit::SegmentLineFitter::HitVec_t = std::vector<Hit_t>

Collection of space points.

Definition at line 37 of file SegmentLineFitter.h.

◆ LinePar_t

using MuonR4::SegmentFit::SegmentLineFitter::LinePar_t = Fitter_t::ParamVec_t

Abrivation of the fitted line parameters.

Definition at line 33 of file SegmentLineFitter.h.

◆ Result_t

using MuonR4::SegmentFit::SegmentLineFitter::Result_t = Fitter_t::FitResult<HitVec_t>

Abrivation of the fit result.

Definition at line 41 of file SegmentLineFitter.h.

◆ Selector_t

using MuonR4::SegmentFit::SegmentLineFitter::Selector_t = Fitter_t::Selector_t<CalibratedSpacePoint>

Abrivation of the hit selector to choose valid hits.

Definition at line 45 of file SegmentLineFitter.h.

Constructor & Destructor Documentation

◆ SegmentLineFitter()

MuonR4::SegmentFit::SegmentLineFitter::SegmentLineFitter	(	const std::string &	name,
		Config &&	config )

Standard constructor.

Parameters

name	Name to be printed in the messaging
config	Fit configuration parameters

Definition at line 79 of file SegmentLineFitter.cxx.

                                                                              :
        AthMessaging{name},
        m_fitter{config, makeActsAthenaLogger(this, name)},
        m_cfg{config} {
        m_goodHitSel.connect<isGoodHit>();
    }

Member Function Documentation

◆ betterResult()

bool MuonR4::SegmentFit::SegmentLineFitter::betterResult	(	const Result_t &	newResult,
		const Result_t &	oldResult ) const

inlineprivate

Returns whether the new fit result is better than the one from the previous iteration.

Selection criterion is the chi2 estimation for the same number of degrees of freedom and then the one which has more degree of freedom but still remains under the good segment threshold

Parameters

newResult	The first fit result
oldResult	The second fit result

Definition at line 417 of file SegmentLineFitter.cxx.

                                                                                                          {
        if (!newResult.converged) {
            ATH_MSG_VERBOSE(__func__<<"() "<<__LINE__<<" The new result did not converge");
            return false;
        }
        const double redChi2New = calcRedChi2(newResult);
        const double redChi2Old = calcRedChi2(oldResult);
        ATH_MSG_VERBOSE(__func__<<"() "<<__LINE__<<" Compare results -- old chi2: "<<redChi2Old<<", nDoF: "
                    <<oldResult.nDoF<<" vs. new chi2: "<<redChi2New<<", nDoF: "<<newResult.nDoF
                    <<" -- outlier removal: "<<m_cfg.outlierRemovalCut);
        if (newResult.nDoF == oldResult.nDoF) {
            return redChi2New < redChi2Old;
        }
        return (redChi2New < m_cfg.outlierRemovalCut && newResult.nDoF > oldResult.nDoF) ||
               (redChi2New > m_cfg.outlierRemovalCut && redChi2New < redChi2Old);
    }

◆ callLineFit()

Result_t MuonR4::SegmentFit::SegmentLineFitter::callLineFit	(	const Acts::CalibrationContext &	cctx,
		const Parameters &	startPars,
		const Amg::Transform3D &	localToGlobal,
		HitVec_t &&	calibHits ) const

private

Calls the underlying line fitter to determine the segment parameters.

Parameters

cctx	Calibration context to fetch later the measurement's calib constants from StoreGate (It's a packed EventContext*)
startPars	Initial line parameters guess
localToGlobal	Transform to align the segment's station inside ATLAS (Mainly neede if the time is fit)
calibHits	List of hits that will be fitted

Check whether a beamspot constraint should be appended

placeholder for a very generous beam spot: 300mm in X,Y (tracking volume), 20000 along Z

Recall that the time is not the same in Acts & Athena

Fit the measurements

Convert back to athena time units

Cache the chi2 terms of the measurements w.r.t. the segment

Definition at line 85 of file SegmentLineFitter.cxx.

                                                                        {
        bool appendsBS = m_cfg.doBeamSpot && countPhiHits(calibHits) > 0;
        
 
        Result_t result{};
        //check the degrees of freedom before try the fit
        if (const std::size_t nPars = m_fitter.config().parsToUse.size(); nPars > 0ul) { 
            auto dOF = m_fitter.countDoF(calibHits, m_goodHitSel);     
            if (dOF.bending + dOF.nonBending < nPars) {
                return result;
            }
            // check that there are at least two crossing stereo measurements
            if (dOF.nonBending == 0ul && nPars == 4ul){
                bool foundU{false}, foundV{false};
                for (const HitVec_t::value_type& hit : calibHits) {
                    if (hit->type() != xAOD::UncalibMeasType::MMClusterType || !isGoodHit(*hit)) {
                        continue;
                    }
                    const auto* mmClust = dynamic_cast<const xAOD::MMCluster*>(hit->spacePoint()->primaryMeasurement());
                    assert(mmClust != nullptr);
                    const auto& design = mmClust->readoutElement()->stripLayer(mmClust->layerHash()).design();
                    if (!design.hasStereoAngle()) {
                        continue;
                    }
                    if (design.stereoAngle() > 0.) {
                        foundU = true;
                    } else {
                        foundV = true;
                    }
                    if (foundU && foundV) {
                        break;
                    }
                }
                if (!foundU || !foundV) {
                    result.measurements = std::move(calibHits);
                    result.parameters = startPars;
                    return result;
                }
                if (m_cfg.doBeamSpot) {
                    appendsBS = true;
                }
            } 
        }
        if (appendsBS) {
            const Amg::Transform3D globToLoc{localToGlobal.inverse()};
            Amg::Vector3D beamSpot{globToLoc.translation()};
            Amg::Vector3D beamLine{globToLoc.linear().col(2)};
            SpacePoint::Cov_t covariance{};
            covariance[toUnderlying(AxisDefs::etaCov)] = square(m_cfg.beamSpotRadius);
            covariance[toUnderlying(AxisDefs::phiCov)] = square(m_cfg.beamSpotLength);
            auto beamSpotSP = std::make_unique<CalibratedSpacePoint>(nullptr, std::move(beamSpot));
            beamSpotSP->setBeamDirection(std::move(beamLine));
            beamSpotSP->setCovariance(std::move(covariance));
            ATH_MSG_VERBOSE(__func__<<"() - "<<__LINE__<<": Beam spot constraint "
                            <<Amg::toString(beamSpotSP->localPosition())<<", "<<beamSpotSP->covariance());
            calibHits.emplace_back(std::move(beamSpotSP));
        }
        ATH_MSG_VERBOSE(__func__<<"() - "<<__LINE__ <<": Start segment fit with parameters "
                <<toString(startPars) <<", plane location: "<<Amg::toString(localToGlobal)
                <<std::endl<<print(calibHits));
 
        FitOpts_t fitOpts{};
        fitOpts.calibContext = cctx;
        fitOpts.calibrator = m_cfg.calibrator;
        fitOpts.selector = m_goodHitSel;
 
        fitOpts.measurements = std::move(calibHits);
        fitOpts.localToGlobal = localToGlobal;
        fitOpts.startParameters = startPars;
        constexpr auto t0idx = toUnderlying(ParamDefs::t0);
        fitOpts.startParameters[t0idx] = ActsTrk::timeToActs(fitOpts.startParameters[t0idx]);
        result = m_fitter.fit(std::move(fitOpts));
        if (m_fitter.config().fitT0) {
            result.parameters[t0idx] = ActsTrk::timeToAthena(result.parameters[t0idx]);
            result.covariance(t0idx, t0idx) = Acts::square(ActsTrk::timeToAthena(1.)) * result.covariance(t0idx, t0idx);
            for (ParamDefs p : {ParamDefs::x0, ParamDefs::y0, ParamDefs::phi, ParamDefs::theta}) {
                auto pidx = toUnderlying(p);
                result.covariance(t0idx, pidx) = ActsTrk::timeToAthena(result.covariance(t0idx, pidx));
                result.covariance(pidx, t0idx) = ActsTrk::timeToAthena(result.covariance(pidx, t0idx));
            }
        }
        {
            const auto[segPos, segDir] = makeLine(result.parameters);
            for (Hit_t& hit : result.measurements) {
                hit->setChi2Term(SeedingAux::chi2Term(segPos, segDir, *hit));
            }
        }
        return result;
    }

◆ checkPrecHitCount()

bool MuonR4::SegmentFit::SegmentLineFitter::checkPrecHitCount ( const HitVec_t & candidateHits ) const

inlineprivate

Checks if the candidate has enough precision hits to fit a segment.

In case of NSW, we check also the orientation of the strips to ensure they provide indipendent constraints on the segment.

Parameters

candidateHits	List of hits on the candidate segment
seed	Segment seed from which the segment was built. Needed for station index

Check whether there is at least one of each micromega strip type. To have a sane topology we need to have at least 2 strips from one kind.

Definition at line 590 of file SegmentLineFitter.cxx.

                                                                                        {
        using namespace Muon::MuonStationIndex;
 
        const size_t nPrecHits = countPrecHits(candidateHits);
        if (nPrecHits < m_cfg.nPrecHitCut) {
            ATH_MSG_VERBOSE(__func__<<"() - "<<__LINE__<<": Not enough precision hits for segment fit. "
                <<nPrecHits<<" < "<<m_cfg.nPrecHitCut);
            return false;
        }
 
        const auto firstHit {std::ranges::find_if(candidateHits, [](const Hit_t& hit){
            return hit->spacePoint() != nullptr;
        })};
        assert(firstHit != candidateHits.end());
        if (toStationIndex((*firstHit)->spacePoint()->msSector()->chamberIndex()) == StIndex::EI &&
            std::ranges::any_of(candidateHits, [](const Hit_t& hit){
                return xAOD::isNSW(hit->type()); })) {
 
            std::array<std::size_t, 3> nStrips{Acts::filledArray<std::size_t, 3>(0u)};
            std::size_t nPhiHits {0u};
            for (const Hit_t& hit : candidateHits) {
                if (!isGoodHit(*hit)) {
                    continue;
                }
                
                if (hit->type() == xAOD::UncalibMeasType::sTgcStripType) {
                    nStrips[0] += isPrecisionHit(*hit);
                    nPhiHits += hit->measuresPhi();
                    continue;
                } else if (hit->type() == xAOD::UncalibMeasType::MMClusterType) {
                    const auto* mmClust = dynamic_cast<const xAOD::MMCluster*>(hit->spacePoint()->primaryMeasurement());
                    assert(mmClust);
                    const auto& design = mmClust->readoutElement()->stripLayer(mmClust->measurementHash()).design();
                    if (!design.hasStereoAngle()) {
                        ++nStrips[0];
                    } else if (design.stereoAngle() > 0.) {
                        ++nStrips[1];
                    } else {
                        ++nStrips[2];
                    }
                }
            }
 
            std::size_t nEtaOrientations = 
                std::ranges::count_if(nStrips, [](std::size_t n){ return n > 0; });
            if (nEtaOrientations == 3u) {
                nEtaOrientations += std::ranges::any_of(  nStrips, [](std::size_t n){ return n > 1; });
            }
            ATH_MSG_VERBOSE(__func__<<"() - "<<__LINE__<<":  nHits: "<<candidateHits.size()
                <<", nPhiHits: "<<nPhiHits<<", nEtaOrientations: "<<nEtaOrientations
                <<", N X-strips: "<<nStrips[0]<<", U-strips: "<<nStrips[1]<<", V-strips: "<<nStrips[2]);
 
            if ( nEtaOrientations == 4u ||
                (nEtaOrientations == 3u && nPhiHits >= 1u) ||
                (nEtaOrientations == 2u && nPhiHits >= 2u)|| 
                (std::ranges::any_of(nStrips, [](std::size_t n){ return n >= 2u; }) && nPhiHits >= 2u)) {
                return true;
            }
            return false;
        }
        return true;
    }

◆ cleanStripLayers()

void MuonR4::SegmentFit::SegmentLineFitter::cleanStripLayers ( HitVec_t & hits ) const

inlineprivate

Marks duplicate hits on a strip layer as outliers to avoid competing contributions from the same layers in the fit.

Hits on the same layer are sorted by their chi2 and the worse ones are rejected if they don't provide additional information

Parameters

hits	List of hit measurements to clean

We need to sort out strip hits on the same layer

Loop over the hits to mark the less compatible hits on the layer as outlier

Both hits measure eta. They've been sorted by lower chi2 -> reject b

Definition at line 349 of file SegmentLineFitter.cxx.

                                                                        {
        const SpacePointPerLayerSorter sorter{};
        std::ranges::sort(hits, [&](const Hit_t&a ,const Hit_t& b){
            if (a->isStraw() || b->isStraw()) {
                return !a->isStraw();
            }
            if (a->type() == xAOD::UncalibMeasType::Other || 
                b->type() == xAOD::UncalibMeasType::Other) {
                return a->type() != xAOD::UncalibMeasType::Other;
            }
            const unsigned lay_a = sorter.sectorLayerNum(*a->spacePoint());
            const unsigned lay_b = sorter.sectorLayerNum(*b->spacePoint());
            if (lay_a != lay_b) {
                return lay_a < lay_b;
            }
            const double chi2a = a->chi2Term();
            const double chi2b = b->chi2Term();
            /* Do not accept pad hits even though they've smaller chi2
             * than the neighbouring strip */
            if (a->type() == xAOD::UncalibMeasType::sTgcStripType) {
                const auto* sTgcA = static_cast<const xAOD::sTgcMeasurement*>(a->spacePoint()->primaryMeasurement());
                const auto* sTgcB = static_cast<const xAOD::sTgcMeasurement*>(b->spacePoint()->primaryMeasurement());
                if (sTgcA->channelType() == xAOD::sTgcMeasurement::sTgcChannelTypes::Pad &&
                    sTgcB->channelType() == xAOD::sTgcMeasurement::sTgcChannelTypes::Strip) {
                    return chi2b > m_cfg.recoveryPull;
                } else if (sTgcB->channelType() == xAOD::sTgcMeasurement::sTgcChannelTypes::Pad &&
                           sTgcA->channelType() == xAOD::sTgcMeasurement::sTgcChannelTypes::Strip) {
                    return chi2a < m_cfg.recoveryPull;
                }
            }
            return chi2a < chi2b;
        });
 
        ATH_MSG_VERBOSE(__func__<<"() - "<<__LINE__ <<": Check for duplicate strip hits");
        for (HitVec_t::iterator itr = hits.begin(); itr != hits.end(); ++itr) {
            const Hit_t& hit_a{*itr};
            if (hit_a->isStraw()){
                break;
            }
            if(hit_a->fitState() == HitState::Duplicate || 
               hit_a->type() == xAOD::UncalibMeasType::Other) {
                continue;
            }
            const unsigned lay_a = sorter.sectorLayerNum(*hit_a->spacePoint());
            for (HitVec_t::iterator itr2 = itr + 1; itr2 != hits.end(); ++itr2) {
                const Hit_t& hit_b{*itr2};
                if (hit_b->type() == xAOD::UncalibMeasType::Other ||
                    hit_b->fitState() == HitState::Duplicate) {
                    continue;
                }
                if (lay_a != sorter.sectorLayerNum(*hit_b->spacePoint())) {
                    break;
                }
                if ((hit_a->measuresEta() && hit_b->measuresEta()) ||
                    (hit_a->measuresPhi() && hit_b->measuresPhi())) {
                    ATH_MSG_VERBOSE(__func__<<"() - "<<__LINE__ <<": Duplicate hit on same layer"<<std::endl
                        <<" -- reject: "<<(*hit_b)<<std::endl
                        <<" -- accept: "<<(*hit_a));
                    hit_b->setFitState(HitState::Duplicate);
                }
            }
        }
    }

◆ convertToSegment()

std::unique_ptr< Segment > MuonR4::SegmentFit::SegmentLineFitter::convertToSegment	(	const Amg::Transform3D &	locToGlobTrf,
		const SegmentSeed *	parentSeed,
		Result_t &&	toConvert ) const

private

Converts the fit result into a segment object.

Parameters

locToGlobTrf	Local to global transform to translate the segment parameters into global parameters
parentSeed	Segment seed from which the segment was built
toConvert	Fitted segment that needs conversion

Definition at line 223 of file SegmentLineFitter.cxx.

                                                                   {
        const auto [locPos, locDir] = makeLine(data.parameters);
        Amg::Vector3D globPos = locToGlob * locPos;
        Amg::Vector3D globDir = locToGlob.linear()* locDir;
 
        std::ranges::sort(data.measurements, [](const Hit_t& a, const Hit_t& b){
            return a->localPosition().z() < b->localPosition().z(); 
        });
        ATH_MSG_VERBOSE(__func__<<"() - "<<__LINE__ <<": Create new segment "
                        <<toString(data.parameters)<<" in "<<patternSeed->msSector()->identString()
                        <<"built from:\n"<<print(data.measurements));
 
        auto finalSeg = std::make_unique<Segment>(std::move(globPos), std::move(globDir),
                                                  patternSeed, std::move(data.measurements),
                                                  data.chi2, data.nDoF);
        finalSeg->setCallsToConverge(data.nIter);
        finalSeg->setParUncertainties(std::move(data.covariance));
        if (m_fitter.config().fitT0) {
            finalSeg->setSegmentT0(data.parameters[toUnderlying(ParamDefs::t0)]);
        }
        return finalSeg;
    }

◆ eraseWrongHits()

void MuonR4::SegmentFit::SegmentLineFitter::eraseWrongHits ( Result_t & candidate ) const

private

Removes all hits from the segment which are obvious outliers.

E.g. tubes which cannot be crossed by the segment.

Parameters

candidate Reference of the segment candidate to prune.

The segment has never crossed the tube

Definition at line 327 of file SegmentLineFitter.cxx.

                                                                    {
        auto [segPos, segDir] = makeLine(candidate.parameters); 
        cleanStripLayers(candidate.measurements);
        candidate.measurements.erase(std::remove_if(candidate.measurements.begin(), 
                                                    candidate.measurements.end(),
            [&](const HitVec_t::value_type& hit){
                if (hit->fitState() == HitState::Valid) {
                    return false;
                } else if (hit->fitState() == HitState::Duplicate) {
                    return true;
                }
                if (hit->type() == xAOD::UncalibMeasType::MdtDriftCircleType) {
                    const double dist = Amg::lineDistance(segPos, segDir, 
                                                            hit->localPosition(), 
                                                            hit->sensorDirection());
                    const auto* dc = static_cast<const xAOD::MdtDriftCircle*>(hit->spacePoint()->primaryMeasurement());
                    return dist >= dc->readoutElement()->innerTubeRadius();
                }
                return false;
            }), candidate.measurements.end());
    }

◆ fitSegment()

std::unique_ptr< Segment > MuonR4::SegmentFit::SegmentLineFitter::fitSegment	(	const EventContext &	ctx,
		const SegmentSeed *	parent,
		const LinePar_t &	startPars,
		const Amg::Transform3D &	localToGlobal,
		HitVec_t &&	calibHits ) const

Fit a set of measurements to a straight segment line.

Badish initial fits are cleaned and then holes are put filled back Returns a nullptr if the fit failed

Parameters

ctx	EventContext to access the calibration constants
parent	Pointer to the seed from which the hits to fit are taken. The seed gives also access to the parent bucket to recover lost hits
startPars	List of parameters serving as an initial guess
localToGlobal	Transform to align the segment's station inside ATLAS (Mainly neede if the time is fit)
calibHits	List of hits that will be fitted

Definition at line 185 of file SegmentLineFitter.cxx.

                                                                  {
        
        const Acts::CalibrationContext cctx = ActsTrk::getCalibrationContext(ctx);
        if (!checkPrecHitCount(calibHits) ) {
            ATH_MSG_VERBOSE(__func__<<"() - "<<__LINE__ <<": Not enough degree of freedom available. What shall be fitted?!");
            return nullptr;
        }
        if (m_cfg.visionTool) {
            Result_t preFit{};
            preFit.parameters = startPars;
            preFit.measurements = calibHits;
            auto seedCopy = convertToSegment(localToGlobal, parent, std::move(preFit));
            m_cfg.visionTool->visualizeSegment(ctx, *seedCopy, "Prefit");
        }
        Result_t segFit = callLineFit(cctx, startPars, localToGlobal, std::move(calibHits));
        if (m_cfg.visionTool && segFit.converged) {
            auto seedCopy = convertToSegment(localToGlobal, parent, Result_t{segFit});
            m_cfg.visionTool->visualizeSegment(ctx, *seedCopy, "Intermediate fit"); 
        }
        if (!removeOutliers(cctx, *parent, localToGlobal,
                            segFit.converged? segFit.parameters : startPars,
                            segFit)) {
            return nullptr;
        }          
        if (!plugHoles(cctx, *parent, localToGlobal, segFit)) {           
            return nullptr;
        }
        auto finalSeg = convertToSegment(localToGlobal, parent, std::move(segFit));
        if (m_cfg.visionTool) {
            m_cfg.visionTool->visualizeSegment(ctx, *finalSeg, "Final fit");
        }
        return finalSeg;
    }

◆ initMessaging()

void AthMessaging::initMessaging ( ) const

privateinherited

Initialize our message level and MessageSvc.

This method should only be called once.

Definition at line 39 of file AthMessaging.cxx.

{
  m_imsg = Athena::getMessageSvc();
  // If user did not set an explicit level, set a default
  if (m_lvl == MSG::NIL) {
    m_lvl = m_imsg ?
      static_cast<MSG::Level>( m_imsg.load()->outputLevel(m_nm) ) :
      MSG::INFO;
  }
}

◆ msg() [1/2]

MsgStream & AthMessaging::msg ( ) const

inlineinherited

The standard message stream.

Returns a reference to the default message stream May not be invoked before sysInitialize() has been invoked.

Definition at line 167 of file AthMessaging.h.

{
  MsgStream* ms = m_msg_tls.get();
  if (!ms) {
    if (!m_initialized.test_and_set()) initMessaging();
    ms = new MsgStream(m_imsg,m_nm);
    m_msg_tls.reset( ms );
  }
 
  ms->setLevel (m_lvl);
  return *ms;
}

◆ msg() [2/2]

MsgStream & AthMessaging::msg ( const MSG::Level lvl ) const

inlineinherited

The standard message stream.

Returns a reference to the default message stream May not be invoked before sysInitialize() has been invoked.

Definition at line 182 of file AthMessaging.h.

183{ return msg() << lvl; }

AthMessaging::msg

MsgStream & msg() const

The standard message stream.

Definition AthMessaging.h:167

◆ msgLvl()

bool AthMessaging::msgLvl ( const MSG::Level lvl ) const

inlineinherited

Test the output level.

Parameters

lvl	The message level to test against

Returns: boolean Indicating if messages at given level will be printed

Return values

true	Messages at level "lvl" will be printed

Definition at line 151 of file AthMessaging.h.

{
  // If user did not set explicit message level we have to initialize
  // the messaging and retrieve the default via the MessageSvc.
  if (m_lvl==MSG::NIL && !m_initialized.test_and_set()) initMessaging();
 
  if (m_lvl <= lvl) {
    msg() << lvl;
    return true;
  } else {
    return false;
  }
}

◆ plugHoles()

bool MuonR4::SegmentFit::SegmentLineFitter::plugHoles	(	const Acts::CalibrationContext &	cctx,
		const SegmentSeed &	seed,
		const Amg::Transform3D &	localToGlobal,
		Result_t &	toRecover ) const

private

Recovery of missed hits.

Hits in the space point bucket that are maximally <RecoveryPull> away from the fitted segment are put onto the segment candidate and the candidate is refitted. If the refitted candidate has a chi2/nDoF < <OutlierRemoval> the canidate is automatically choosen otherwise, its chi needs to be better.

Parameters

cctx	Calibration context to fetch later the measurement's calib constants from StoreGate (It's a packed EventContext*)
seed	Parent seed from which the segment fit is actually triggered The seed is mainly used for visualization purposes
localToGlobal	Transform to align the spectrometer sector within ATLAS mainly used for the t0 fit
fitResult	Previously achieved fit result to be checked. The measurements on the result and the paramters are updated accordingly

We've the first estimator of the segment fit

Setup a map to replace space points if they better suite

Loop over all hits in the parent bucket

Hit already used in the segment fit

If the hit is a phi measurement check at least if it can be hit by the segment

Use the pull of the uncalibrated measurement to estimate whether a calibration is actually worth

No extra hit has been found

Remove the beamspot constraint measurement

If the chi2 is less than 5, no outlier rejection is launched. So also accept any recovered segment below that threshold

Next check whether the recovery made measurements marked as outlier feasable to the hole recovery

Definition at line 433 of file SegmentLineFitter.cxx.

                                                                 {
        ATH_MSG_VERBOSE(__func__<<"() - "<<__LINE__ <<": segment "<<toString(toRecover.parameters)
                        <<", chi2: "<< calcRedChi2(toRecover) <<", nDoF: "<<toRecover.nDoF);
        
 
        std::unordered_set<const SpacePoint*> usedSpacePoints{};
        for (auto& hit : toRecover.measurements) {
            ATH_MSG_VERBOSE(__func__<<"() - "<<__LINE__ <<": "<<(*hit)<<" is known");
            usedSpacePoints.insert(hit->spacePoint());
        }
        const EventContext& ctx{*cctx.get<const EventContext*>()};
        
        const double timeOff = toRecover.parameters[toUnderlying(ParamDefs::t0)];
        HitVec_t candidateHits{};
        std::size_t hasCandidate{0};
        const auto [locPos, locDir] = makeLine(toRecover.parameters);
        for (const auto& hit : *seed.parentBucket()){            
            if (usedSpacePoints.count(hit.get())){ 
                continue;
            }
            Hit_t calibHit{};
            double pull{-1.};
            if (hit->isStraw()) {
                using namespace Acts::detail::LineHelper;
                const double dist = signedDistance(locPos, locDir, hit->localPosition(), hit->sensorDirection());
                const auto* dc = static_cast<const xAOD::MdtDriftCircle*>(hit->primaryMeasurement());
                // Check whether the tube is crossed by the hit 
                if (std::abs(dist) >= dc->readoutElement()->innerTubeRadius()) {
                    continue;
                }
            } else {
                if (!hit->measuresEta() && 
                    std::abs(hit->sensorDirection().dot(hit->localPosition() - 
                        SeedingAux::extrapolateToPlane(locPos,locDir, *hit))) >
                        1.1*std::sqrt(hit->covariance()[toUnderlying(AxisDefs::etaCov)])){
                    continue;
                }
                pull = std::sqrt(SeedingAux::chi2Term(locPos, locDir, *hit));
                if (pull > 1.1 * m_cfg.recoveryPull) {
                    continue;
                }
            }
            calibHit = m_cfg.calibrator->calibrate(ctx, hit.get(), locPos, locDir, ActsTrk::timeToActs(timeOff));
            calibHit->setChi2Term(SeedingAux::chi2Term(locPos, locDir, *calibHit));
            if (calibHit->chi2Term() <= Acts::square(m_cfg.recoveryPull)) {
                hasCandidate += calibHit->fitState() == HitState::Valid;
                ATH_MSG_VERBOSE(__func__<<"() - "<<__LINE__<<": Candidate hit for recovery "
                            <<(*calibHit));
            } else {
                calibHit->setFitState(HitState::Outlier);
                ATH_MSG_VERBOSE(__func__<<"() - "<<__LINE__<<": Outlier hit "
                            <<(*calibHit)<<" -> limit: "<<m_cfg.recoveryPull);
            }
            candidateHits.push_back(std::move(calibHit));                
        }
        if (!hasCandidate) {
            ATH_MSG_VERBOSE(__func__<<"() - "<<__LINE__<<": No space point candidates for recovery were found");
            toRecover.measurements.insert(toRecover.measurements.end(), 
                                          std::make_move_iterator(candidateHits.begin()),
                                          std::make_move_iterator(candidateHits.end()));
            eraseWrongHits(toRecover);
            return toRecover.nDoF > 0;
        }
        ATH_MSG_VERBOSE(__func__<<"() - "<<__LINE__<<": Found "<<hasCandidate<<" space points for recovery. ");
 
 
        HitVec_t hitsForRecovery = toRecover.measurements;
        if (m_cfg.doBeamSpot) {
            removeBeamSpot(hitsForRecovery);
        }
 
        hitsForRecovery.insert(hitsForRecovery.end(), 
                               candidateHits.begin(),
                               candidateHits.end());
 
        cleanStripLayers(hitsForRecovery);
 
        Result_t recovered = callLineFit(cctx, toRecover.parameters, localToGlobal, 
                                         std::move(hitsForRecovery));
       
        if (betterResult(recovered, toRecover)) {
            ATH_MSG_VERBOSE(__func__<<"() - "<<__LINE__<<": Accept segment with recovered "
                            <<(recovered.nDoF  - toRecover.nDoF)<<" extra nDoF.");
            recovered.nIter += toRecover.nIter;
            toRecover = std::move(recovered);
 
            std::vector<const CalibratedSpacePoint*> stripOutliers{};
            stripOutliers.reserve(toRecover.measurements.size());
            unsigned recovLoop{(candidateHits.size() != hasCandidate)*m_cfg.nRecoveryLoops};
            while (++recovLoop <= m_cfg.nRecoveryLoops) {   
                ATH_MSG_VERBOSE(__func__<<"() - "<<__LINE__<<": Enter recovery loop "<<recovLoop<<".");
                hitsForRecovery = toRecover.measurements;
                // Remove the beamspot
                if (m_cfg.doBeamSpot) {
                    removeBeamSpot(hitsForRecovery);
                }
                // Check whether an outlier can be lifted to on-track
                for (HitVec_t::value_type& hit : hitsForRecovery) {
                    if (hit->fitState() != HitState::Outlier) {
                        continue;
                    }
                    if (hit->chi2Term() < Acts::square(m_cfg.recoveryPull)) {
                        ATH_MSG_VERBOSE(__func__<<"() - "<<__LINE__<<": Try to recover outlier "<<(*hit));
                        hit->setFitState(HitState::Valid);
                        stripOutliers.push_back(hit.get());
                    } 
                }
                // Nothing to recover
                if (stripOutliers.empty()) {
                    ATH_MSG_VERBOSE(__func__<<"() - "<<__LINE__<<": No additional measurement found");
                    break;
                }
                // Ensure that only one hit per layer is fit
                cleanStripLayers(hitsForRecovery);
                // Recovery turned out to be duplicates on the same layer
                if (std::ranges::none_of(stripOutliers,[](const CalibratedSpacePoint* sp){
                        return sp->fitState() == HitState::Valid;
                    })) {
                    ATH_MSG_VERBOSE(__func__<<"() - "<<__LINE__<<": Outliers turned out to be duplicates.");
                    break;
                }
                ATH_MSG_VERBOSE(__func__<<"() - "<<__LINE__<<": Start fit without the outliers.");
                stripOutliers.clear();
                recovered = callLineFit(cctx, toRecover.parameters, localToGlobal, std::move(hitsForRecovery));
                if (!betterResult(recovered, toRecover)) {
                    break;
                }
                recovered.nIter += toRecover.nIter;
                toRecover = std::move(recovered);
            }
        } else{
            for (HitVec_t::value_type& hit : candidateHits) {
                hit->setFitState(HitState::Outlier);
                toRecover.measurements.push_back(std::move(hit));
            }
        }
        eraseWrongHits(toRecover);
        return true;
    }        

◆ removeOutliers()

bool MuonR4::SegmentFit::SegmentLineFitter::removeOutliers	(	const Acts::CalibrationContext &	cctx,
		const SegmentSeed &	seed,
		const Amg::Transform3D &	localToGlobal,
		const LinePar_t &	startPars,
		Result_t &	fitResult ) const

private

Cleans the fitted segment from the most outlier hit and then attempts to refit the segment.

The outlier removal is not run if the segment has already a chi2 / nDoF better than <outlierRemovalCut>. Returns false if the recovery lead to the destruction of all nDoF

Parameters

cctx	Calibration context to fetch later the measurement's calib constants from StoreGate (It's a packed EventContext*)
seed	Parent seed from which the segment fit is actually triggered The seed is mainly used for visualization purposes
localToGlobal	Transform to align the spectrometer sector within ATLAS mainly used for the t0 fit
startPars	The initial parameters from which the fit shall be launched In case of out of bound parameters, the start parameters are returned otherwise the last obtained fit parameters
fitResult	Previously achieved fit result to be checked. The measurements on the result and the paramters are updated accordingly

Remove a priori the beamspot constaint as it never should pose any problem and another one will be added anyway in the next iteration

Next sort the measurements by ascending chi2

Move the outliers to the front

Check again the available DOF and number of precision hits after hit removal

Refit the segment line without the measurement

Definition at line 248 of file SegmentLineFitter.cxx.

                                                                      {
 
        if (!checkPrecHitCount(fitResult.measurements) || 
            fitResult.nIter > m_fitter.config().maxIter) {
            ATH_MSG_VERBOSE(__func__<<"() - "<<__LINE__ 
                            <<": No degree of freedom available. What shall be removed?!. nDoF: "
                            <<fitResult.nDoF<<", n-meas: "<<countPrecHits(fitResult.measurements)
                            <<std::endl<<print(fitResult.measurements));
            return false;
        }
        if (fitResult.converged && calcRedChi2(fitResult) < m_cfg.outlierRemovalCut) {
            ATH_MSG_VERBOSE(__func__<<"() - "<<__LINE__ <<": The segment "<<toString(fitResult.parameters)
                          <<" is already of good quality "<< calcRedChi2(fitResult)<<". Don't remove outliers");
            return true;
        }
        if (fitResult.nDoF == 0u){
            return false;
        }
        ATH_MSG_VERBOSE(__func__<<"() - "<<__LINE__ <<": Segment "
                       <<toString(fitResult.parameters)<<", nIter: "<<fitResult.nIter
                       <<" is of badish quality. "<<print(fitResult.measurements)
                       <<std::endl<<"Remove worst hit");
        
        if (m_cfg.doBeamSpot) {
            removeBeamSpot(fitResult.measurements);
        }
        std::ranges::sort(fitResult.measurements,
                [](const HitVec_t::value_type& a, const HitVec_t::value_type& b){
                    if (isGoodHit(*a) != isGoodHit(*b)) {
                        return !isGoodHit(*a);
                    }
                    return a->chi2Term() < b->chi2Term();                   
                });
        fitResult.measurements.back()->setFitState(HitState::Outlier);
        ATH_MSG_VERBOSE(__func__<<"() - "<<__LINE__<<" Mark "<<(*fitResult.measurements.back())<<" as outlier");
        if (!checkPrecHitCount(fitResult.measurements)) {
            ATH_MSG_VERBOSE(__func__<<"() - "<<__LINE__ 
                <<": No degree of freedom available after outlier removal. n-meas: "
                <<countPrecHits(fitResult.measurements)<<std::endl<<print(fitResult.measurements));
            return false;
        }
        Result_t newAttempt = callLineFit(cctx, startPars, localToGlobal, 
                                          std::move(fitResult.measurements));
        if (newAttempt.converged) {
            ATH_MSG_VERBOSE(__func__<<"() - "<<__LINE__<<" The outlier removal converged.");
            newAttempt.nIter+=fitResult.nIter;
            fitResult = std::move(newAttempt);
             if (m_cfg.visionTool) {
                const EventContext& ctx{*cctx.get<const EventContext*>()};
                auto seedCopy = convertToSegment(localToGlobal, &seed, Result_t{fitResult});
                m_cfg.visionTool->visualizeSegment(ctx, *seedCopy, "Bad fit recovery");
            }
        } else {
            ATH_MSG_VERBOSE(__func__<<"() - "<<__LINE__
                <<" Outlier removal fit did not converge. Needed iterations: "<<newAttempt.nIter);
            if (newAttempt.nIter == 0ul) {
                return false;
            }
            fitResult.nIter+=newAttempt.nIter;
            fitResult.measurements = std::move(newAttempt.measurements);
        }
        return removeOutliers(cctx, seed, localToGlobal,
                              fitResult.converged ? fitResult.parameters : startPars, 
                              fitResult);
    }

◆ setLevel()

void AthMessaging::setLevel ( MSG::Level lvl )

inherited

Change the current logging level.

Use this rather than msg().setLevel() for proper operation with MT.

Definition at line 28 of file AthMessaging.cxx.

{
  m_lvl = lvl;
}

Member Data Documentation

◆ ATLAS_THREAD_SAFE

std::atomic_flag m_initialized AthMessaging::ATLAS_THREAD_SAFE = ATOMIC_FLAG_INIT

mutableprivateinherited

Messaging initialized (initMessaging).

Definition at line 141 of file AthMessaging.h.

◆ m_cfg

ConfigSwitches MuonR4::SegmentFit::SegmentLineFitter::m_cfg {}

private

Configuration switches of the ATLAS fitter implementation.

Definition at line 108 of file SegmentLineFitter.h.

108{};

◆ m_fitter

Fitter_t MuonR4::SegmentFit::SegmentLineFitter::m_fitter

private

Actual implementation of the straight line fit.

Definition at line 106 of file SegmentLineFitter.h.

◆ m_goodHitSel

Selector_t MuonR4::SegmentFit::SegmentLineFitter::m_goodHitSel {}

private

Selector to identify the valid hits.

Definition at line 110 of file SegmentLineFitter.h.

110{};

◆ m_imsg

std::atomic<IMessageSvc*> AthMessaging::m_imsg { nullptr }

mutableprivateinherited

MessageSvc pointer.

Definition at line 135 of file AthMessaging.h.

135{ nullptr };

◆ m_lvl

std::atomic<MSG::Level> AthMessaging::m_lvl { MSG::NIL }

mutableprivateinherited

Current logging level.

Definition at line 138 of file AthMessaging.h.

138{ MSG::NIL };

◆ m_msg_tls

boost::thread_specific_ptr<MsgStream> AthMessaging::m_msg_tls

mutableprivateinherited

MsgStream instance (a std::cout like with print-out levels).

Definition at line 132 of file AthMessaging.h.

◆ m_nm

std::string AthMessaging::m_nm

privateinherited

Message source name.

Definition at line 129 of file AthMessaging.h.

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

Classes

Public Types

Public Member Functions

Private Member Functions

Private Attributes

Detailed Description

Member Typedef Documentation

◆ FitOpts_t

◆ FitPars_t

◆ Fitter_t

◆ Hit_t

◆ HitState

◆ HitVec_t

◆ LinePar_t

◆ Result_t

◆ Selector_t

Constructor & Destructor Documentation

◆ SegmentLineFitter()

Member Function Documentation

◆ betterResult()

◆ callLineFit()

◆ checkPrecHitCount()

◆ cleanStripLayers()

◆ convertToSegment()

◆ eraseWrongHits()

◆ fitSegment()

◆ initMessaging()

◆ msg() [1/2]

◆ msg() [2/2]

◆ msgLvl()

◆ plugHoles()

◆ removeOutliers()

◆ setLevel()

Member Data Documentation

◆ ATLAS_THREAD_SAFE

◆ m_cfg

◆ m_fitter

◆ m_goodHitSel

◆ m_imsg

◆ m_lvl

◆ m_msg_tls

◆ m_nm