SegmentAmbiSolver.cxx

Go to the documentation of this file.

/*
  Copyright (C) 2002-2024 CERN for the benefit of the ATLAS collaboration
*/
#include <MuonPatternHelpers/SegmentAmbiSolver.h>
#include <MuonPatternHelpers/SegmentFitHelperFunctions.h>
#include <Acts/Utilities/Enumerate.hpp>
 
namespace MuonR4 {
    using namespace 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() / segment.nDoF();
    }
    SegmentVec SegmentAmbiSolver::resolveAmbiguity(const ActsGeometryContext& 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 (std::unique_ptr<Segment>& resolveMe : toResolve) {
 
            MeasurementSet testMeas{extractPrds(*resolveMe)};
 
            Resolution reso{Resolution::noOverlap};
            unsigned int resolvedIdx{0};
            for (std::unique_ptr<Segment>& goodSeg : resolved) {
                ATH_MSG_VERBOSE("Test against segment "<<toString(localSegmentPars(gctx, *goodSeg)));
                MeasurementSet& resolvedM = segMeasurements[resolvedIdx];
                std::vector<int>& existSigns{segmentSigns[resolvedIdx++]};
                unsigned int 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 int sameSides{0};
                for (unsigned int s =0 ; s < existSigns.size(); ++s) {
                    sameSides += (reEvaluatedSigns[s] == existSigns[s]);
                }
                if (sameSides != existSigns.size() && resolveMe->nDoF() == goodSeg->nDoF()) {
                    ATH_MSG_VERBOSE("Reference signs: "<<existSigns<<" / re-evaluated: "<<reEvaluatedSigns);
                    continue;
                }
                reso = Resolution::subSet;
                const double resolvedChi2  = redChi2(*goodSeg);
                const double resolveMeChi2 = redChi2(*resolveMe);
                ATH_MSG_VERBOSE("Chi2 good "<<resolvedChi2<<", candidate chi2: "<<resolveMeChi2);
                if (resolveMeChi2 < m_cfg.selectByNDoFChi2 && resolvedChi2 < m_cfg.selectByNDoFChi2) {
                    if (resolveMe->nDoF() > goodSeg->nDoF()){
                        reso = Resolution::superSet;
                    }
                } else if (resolveMeChi2 < resolvedChi2) {
                    reso = Resolution::superSet;
                }
                if (reso == Resolution::superSet) {
                    std::swap(goodSeg, resolveMe);
                    std::swap(resolvedM, testMeas);
                    existSigns = driftSigns(gctx, *resolveMe, resolveMe->measurements());
                }
            }
            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 ActsGeometryContext& gctx,
                                                   const Segment& segment,
                                                   const Segment::MeasVec& measurements) const {
        std::vector<int> signs{};
        signs.reserve(measurements.size());
        const auto [locPos, locDir] = makeLine(localSegmentPars(gctx, segment));
        ATH_MSG_VERBOSE("Fetch drift signs for segment "<<segment.msSector()->identString()<<" -- "<<Amg::toString(locPos)
                        <<Amg::toString(locDir));
        for (const Segment::MeasVec::value_type& hit : measurements) {
            if (!hit->spacePoint()) {
                continue;
            }
            signs.push_back(SegmentFitHelpers::driftSign(locPos,locDir,*hit, msg()));
        }
        return signs;
    }
    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 int 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);
            }); 
    }
}