SegmentAmbiSolver.cxx

Go to the documentation of this file.

/*
  Copyright (C) 2002-2025 CERN for the benefit of the ATLAS collaboration
*/
#include <MuonPatternHelpers/SegmentAmbiSolver.h>
#include <Acts/Utilities/Enumerate.hpp>
 
namespace MuonR4::SegmentFit {
    using SegmentVec = SegmentAmbiSolver::SegmentVec;
 
    SegmentAmbiSolver::SegmentAmbiSolver(const std::string& name, Config&& cfg):
        AthMessaging{name},
        m_cfg{std::move(cfg)} {}
 
    double SegmentAmbiSolver::redChi2(const Segment& segment) const {
        return segment.chi2() / std::max(segment.nDoF(), 1u);
    }
    SegmentVec SegmentAmbiSolver::resolveAmbiguity(const ActsTrk::GeometryContext& gctx,
                                                   SegmentVec&& toResolve) const {
        
        std::ranges::stable_sort(toResolve,[this](const SegmentVec::value_type& a,
                                                  const SegmentVec::value_type& b){
            const double redChi2A = redChi2(*a);
            const double redChi2B = redChi2(*b);
            if (redChi2A < m_cfg.selectByNDoFChi2 && redChi2B < m_cfg.selectByNDoFChi2){
                return a->nDoF() > b->nDoF();
            }
            return redChi2A < redChi2B;
        });
        SegmentVec resolved{};

        resolved.push_back(std::move(toResolve[0]));
        toResolve.erase(toResolve.begin());
        std::vector<std::vector<int>> segmentSigns{driftSigns(gctx, *resolved.front(), resolved.front()->measurements())};
        std::vector<MeasurementSet> segMeasurements{extractPrds(*resolved.front())};
 
        for (SegmentVec::value_type& resolveMe : toResolve) {
            ATH_MSG_VERBOSE("Try to resolve new segment "<<toString(localSegmentPars(gctx, *resolveMe))
                          <<" redChi2: "<<redChi2(*resolveMe)<<" nDoF: "<<resolveMe->nDoF());
            MeasurementSet testMeas{extractPrds(*resolveMe)};
            Resolution reso{Resolution::noOverlap};
            unsigned resolvedIdx{0};
            for (SegmentVec::value_type& goodSeg : resolved) {
                ATH_MSG_VERBOSE("Test against segment "<<toString(localSegmentPars(gctx, *goodSeg))
                        <<" redChi2: "<<redChi2(*goodSeg)<<" nDoF: "<<goodSeg->nDoF());
                MeasurementSet& resolvedM = segMeasurements[resolvedIdx];
                std::vector<int>& existSigns{segmentSigns[resolvedIdx++]};
                unsigned shared = countShared(resolvedM, testMeas);
                if (shared < m_cfg.sharedPrecHits) {
                    ATH_MSG_VERBOSE("Too few shared measurements "<<shared<<" (Required: "<<m_cfg.sharedPrecHits<<").");
                    continue;
                }
                const std::vector<int> reEvaluatedSigns{driftSigns(gctx, *resolveMe, goodSeg->measurements())};
 
                unsigned sameSides{0};
                for (unsigned s =0 ; s < existSigns.size(); ++s) {
                    sameSides += (reEvaluatedSigns[s] == existSigns[s]);
                }
                if (!m_cfg.remLeftRightAmbi && sameSides != existSigns.size() && resolveMe->nDoF() == goodSeg->nDoF()) {
                    ATH_MSG_VERBOSE("Reference signs: "<<existSigns<<" / re-evaluated: "<<reEvaluatedSigns);
                    continue;
                }
                const double resolvedChi2  = redChi2(*goodSeg);
                const double resolveMeChi2 = redChi2(*resolveMe);
                reso = resolveMeChi2 < resolvedChi2 ? Resolution::superSet : Resolution::subSet;
 
                ATH_MSG_VERBOSE("Chi2 good "<<resolvedChi2<<", candidate chi2: "<<resolveMeChi2);
                if (resolveMeChi2 < m_cfg.selectByNDoFChi2 && resolvedChi2 < m_cfg.selectByNDoFChi2) {
                    reso =  goodSeg->nDoF() > resolveMe->nDoF() ? Resolution::subSet : Resolution::superSet;
                }
                if (reso == Resolution::superSet) {
                    std::swap(goodSeg, resolveMe);
                    std::swap(resolvedM, testMeas);
                    existSigns = driftSigns(gctx, *resolveMe, resolveMe->measurements());
                    break;
                } else if (reso == Resolution::subSet) {
                    break;
                }
            }
            if (reso == Resolution::noOverlap) {
                segMeasurements.push_back(std::move(testMeas));
                segmentSigns.push_back(driftSigns(gctx, *resolveMe, resolveMe->measurements()));
                resolved.push_back(std::move(resolveMe));
            }
        }
        return resolved;
    }
    std::vector<int> SegmentAmbiSolver::driftSigns(const ActsTrk::GeometryContext& gctx,
                                                   const Segment& segment,
                                                   const Segment::MeasVec& measurements) const {
        const auto [pos, dir] = makeLine(localSegmentPars(gctx, segment));
     
        ATH_MSG_VERBOSE("Fetch drift signs for segment "<<segment.msSector()->identString()
                      <<" -- "<<Amg::toString(pos)<<" + " <<Amg::toString(dir));
        return SeedingAux::strawSigns(pos, dir, measurements);
    }
    SegmentAmbiSolver::MeasurementSet 
        SegmentAmbiSolver::extractPrds(const Segment& segment) const {
        
        MeasurementSet meas{};
        for (const Segment::MeasVec::value_type& hit : segment.measurements()) {
            if (!hit->spacePoint() || hit->fitState() != CalibratedSpacePoint::State::Valid || !hit->measuresEta()) {
                continue;
            }
            meas.insert(hit->spacePoint()->primaryMeasurement());
            if (hit->spacePoint()->secondaryMeasurement()) {
                 meas.insert(hit->spacePoint()->secondaryMeasurement());
            }
        }
        return meas;
    }
    unsigned SegmentAmbiSolver::countShared(const MeasurementSet& measSet1, 
                                            const MeasurementSet& measSet2) const {
        if (measSet1.size() > measSet2.size()) {
            return std::count_if(measSet2.begin(),measSet2.end(),[&measSet1](const xAOD::UncalibratedMeasurement* meas){
                return measSet1.count(meas);
            });
        }
        return std::count_if(measSet1.begin(),measSet1.end(),[&measSet2](const xAOD::UncalibratedMeasurement* meas){
                return measSet2.count(meas);
            }); 
    }
}