SegmentFittingAlg.cxx

Go to the documentation of this file.

/*
  Copyright (C) 2002-2026 CERN for the benefit of the ATLAS collaboration
*/
#include "SegmentFittingAlg.h"
 
#include "GeoPrimitives/GeoPrimitivesHelpers.h"
#include "EventPrimitives/EventPrimitivesHelpers.h"
 
#include "MuonPatternHelpers/MdtSegmentSeedGenerator.h"
#include "MuonSpacePoint/SpacePointPerLayerSplitter.h"
#include "MuonTruthHelpers/MuonSimHitHelpers.h"
#include "ActsCalibBase/CalibrationContext.h"
#include "ActsInterop/Logger.h"
 
#include "MuonVisualizationHelpersR4/VisualizationHelpers.h"
 
#include "Acts/Definitions/Units.hpp"
 
#include <format>
 
using namespace Acts;
using namespace Acts::UnitLiterals;
 
namespace MuonR4 {
    using namespace SegmentFit;
    using namespace MuonValR4;
 
    using PrimitiveVec = MuonValR4::IPatternVisualizationTool::PrimitiveVec;
 
    SegmentFittingAlg::~SegmentFittingAlg() = default;
    StatusCode SegmentFittingAlg::initialize() {
        ATH_CHECK(m_geoCtxKey.initialize());
        ATH_CHECK(m_seedKey.initialize());
        ATH_CHECK(m_outSegments.initialize());  
        ATH_CHECK(m_calibTool.retrieve());
        ATH_CHECK(m_idHelperSvc.retrieve());
        ATH_CHECK(m_visionTool.retrieve(EnableTool{!m_visionTool.empty()}));
        SegmentAmbiSolver::Config cfg{};
        m_ambiSolver = std::make_unique<SegmentAmbiSolver>(name(), std::move(cfg));
 
        SegmentLineFitter::Config fitCfg{};
        fitCfg.calibrator = m_calibTool.get();
        fitCfg.visionTool = m_visionTool.get();
        fitCfg.idHelperSvc = m_idHelperSvc.get();
        fitCfg.fitT0 = m_doT0Fit;
        fitCfg.recalibrate = m_recalibInFit;
        fitCfg.useFastFitter = m_useFastFitter;
        fitCfg.fastPreFitter = m_fastPreFitter;
        fitCfg.ignoreFailedPreFit = m_ignoreFailedPreFit;
        fitCfg.useHessian = m_hessianResidual;
 
        fitCfg.doBeamSpot = m_doBeamspotConstraint;
        fitCfg.beamSpotRadius = m_beamSpotR;
        fitCfg.beamSpotLength = m_beamSpotL;
 
        if (m_doBeamspotConstraint) {
            for (const auto* prop :{&m_seedMaxBsR, &m_seedMaxBsL}) {
                if (prop->value().size() != s_stIdxMax || 
                    std::ranges::min(prop->value()) < 0.) {
                    ATH_MSG_ERROR("Invalid configuratnon of "<<(prop));
                    return StatusCode::FAILURE;
                }
                ATH_MSG_DEBUG("Configured "<<(*prop));
            }
        }
        fitCfg.outlierRemovalCut = m_outlierRemovalCut;
        fitCfg.recoveryPull = m_recoveryPull;
        fitCfg.nPrecHitCut = m_precHitCut;
        fitCfg.maxIter = m_maxIter;
        ATH_MSG_DEBUG("Fitter configuration: \n -- "<<m_doT0Fit<<"\n -- "<<m_recalibInFit<<"\n -- "<<m_useFastFitter
                <<"\n -- "<<m_fastPreFitter<<"\n -- "<<m_ignoreFailedPreFit<<"\n -- "<<m_hessianResidual
                <<"\n -- "<<m_doBeamspotConstraint<<"\n -- "<<m_beamSpotR<<"\n -- "<<m_beamSpotL<<"\n -- "<<m_outlierRemovalCut
                <<"\n -- "<<m_recoveryPull<<"\n -- "<<m_precHitCut<<"\n -- "<<m_maxIter);
 
        m_fitter = std::make_unique<SegmentFit::SegmentLineFitter>(name(), std::move(fitCfg));
 
        MdtSegmentSeedGenerator::Config genCfg{};
        genCfg.hitPullCut = m_seedHitChi2;
        genCfg.busyLayerLimit = m_busyLayerLimit;
        genCfg.startWithPattern = m_tryPatternPars;
        ATH_MSG_DEBUG("Seeder configuration: \n - "<<m_tryPatternPars
                    <<"\n - "<<m_seedHitChi2
                    <<"\n - "<<m_busyLayerLimit);  
        m_seeder = std::make_unique<SegmentFit::MdtSegmentSeedGenerator>(genCfg, makeActsAthenaLogger(this, name()));
        genCfg.busyLayerLimit += 2;
        m_seederBEE = std::make_unique<SegmentFit::MdtSegmentSeedGenerator>(genCfg, makeActsAthenaLogger(this, name()));
     
        return StatusCode::SUCCESS;
    }
    StatusCode SegmentFittingAlg::execute(const EventContext& ctx) const {
        const ActsTrk::GeometryContext* gctx{nullptr};
        ATH_CHECK(SG::get(gctx, m_geoCtxKey, ctx));
        const SegmentSeedContainer* segmentSeeds=nullptr; 
        ATH_CHECK(SG::get(segmentSeeds, m_seedKey, ctx));
    
        SG::WriteHandle writeSegments{m_outSegments, ctx};
        ATH_CHECK(writeSegments.record(std::make_unique<SegmentContainer>()));
        SegmentVec_t allSegments{};
        ATH_MSG_VERBOSE("execute() - Start processing " << segmentSeeds->size() << " pattern seeds.");
        for (const SegmentSeed* seed : *segmentSeeds) {
            SegmentVec_t segments = fitSegmentSeed(ctx, *gctx, seed);
             if (m_visionTool.isEnabled() && segments.size() > 1) {
                auto drawFinalReco = [this, &segments, &gctx, &ctx,&seed](const std::string& nameTag) {
                    PrimitiveVec segmentLines{};
                    double yLegend{0.85};
                    segmentLines.push_back(drawLabel(std::format("# segments: {:d}",  segments.size()), 0.2, yLegend, 14));
                    yLegend-=0.04;
                    for (const std::unique_ptr<Segment>& seg : segments) {
                        const Parameters pars = localSegmentPars(*gctx, *seg);
                        const auto [pos, dir] = makeLine(pars);
                        segmentLines.emplace_back(drawLine(pars, -Gaudi::Units::m, Gaudi::Units::m, kRed));
                        std::stringstream signStream{};
                        const auto& cov = seg->covariance();
                        signStream<<std::format("#chi^{{2}}/nDoF: {:.2f} ({:}), ", seg->chi2() / seg->nDoF(), seg->nDoF());
                        signStream<<std::format("y_{{0}}={:.2f}#pm{:.2f}",
                                                pars[toUnderlying(ParamDefs::y0)],
                                                Amg::error(cov, toUnderlying(ParamDefs::y0)))<<", ";
                        signStream<<std::format("#theta={:.2f}#pm{:.2f}^{{#circ}}", 
                                                pars[toUnderlying(ParamDefs::theta)]/ 1._degree,
                                                Amg::error(cov, toUnderlying(ParamDefs::theta)) / 1._degree )<<", ";
                        for (const Segment::MeasType& m : seg->measurements()) {
                            if (m->type() == xAOD::UncalibMeasType::MdtDriftCircleType && 
                                m->fitState() == CalibratedSpacePoint::State::Valid) {
                                signStream<<(SeedingAux::strawSign(pos, dir, *m) == -1 ? "L" : "R");
                            }
                        }
                        segmentLines.push_back(drawLabel(signStream.str(), 0.2, yLegend, 13));
                        yLegend-=0.03;
                    }
 
                    m_visionTool->visualizeBucket(ctx, *seed->parentBucket(), nameTag, std::move(segmentLines));
                };
                drawFinalReco("all segments");
                const unsigned int nBeforeAmbi = segments.size();
                segments = m_ambiSolver->resolveAmbiguity(*gctx, std::move(segments));
                if (nBeforeAmbi != segments.size()) {
                    drawFinalReco("post ambiguity");
                }
            } else if (m_visionTool.isEnabled() && segments.empty() &&
                      std::ranges::any_of(seed->getHitsInMax(),[this](const SpacePoint* hit){
                            return  m_visionTool->isLabeled(*hit);
                      })) {
                m_visionTool->visualizeSeed(ctx, *seed, "Failed fit");
            }
            allSegments.insert(allSegments.end(), std::make_move_iterator(segments.begin()),
                                                  std::make_move_iterator(segments.end()));
        }
        resolveAmbiguities(*gctx, allSegments);
        writeSegments->insert(writeSegments->end(),
                                  std::make_move_iterator(allSegments.begin()),
                                  std::make_move_iterator(allSegments.end()));
        ATH_MSG_VERBOSE("Found in total "<<writeSegments->size()<<" segments. ");
        return StatusCode::SUCCESS; 
    }
 
    SegmentFittingAlg::SegmentVec_t
         SegmentFittingAlg::fitSegmentSeed(const EventContext& ctx,
                                           const ActsTrk::GeometryContext& gctx,
                                           const SegmentSeed* patternSeed) const {
 
        const auto& locToGlob = patternSeed->msSector()->localToGlobalTransform(gctx);
        std::vector<std::unique_ptr<Segment>> segments{};
 
        const Acts::CalibrationContext cctx = ActsTrk::getCalibrationContext(ctx);
 
        using namespace Muon::MuonStationIndex;
        using State_t = MdtSegmentSeedGenerator::State_t;
        // Make sure patternSeed->parameters() are in ACTS units!
        State_t seedState{patternSeed->parameters(), patternSeed, m_calibTool.get(), m_recalibSeed,
                          [&](const Amg::Vector3D& tangentSeedPos, const Amg::Vector3D& tangentSeedDir){
                                if (!m_doBeamspotConstraint) {
                                    return true;
                                }
                                const Amg::Vector3D globPos{locToGlob*tangentSeedPos};
                                const Amg::Vector3D globDir{locToGlob.linear()*tangentSeedDir};

                                const StIndex stIdx = toStationIndex(patternSeed->msSector()->chamberIndex());
                                switch (stIdx) {
                                    using enum StIndex;
                                    case EM: {
                                        if (std::abs(globPos.eta()) > 2.35 && 
                                            std::abs(globPos.theta() - globDir.theta()) < 5.8_degree) {
                                            return true;
                                        }
                                        break;
                                    }
                                    case EO: {
                                         if (std::abs(globPos.eta()) > 2.35 && 
                                            std::abs(globPos.theta() - globDir.theta()) < 3._degree) {
                                            return true;
                                        }
                                        break;
                                    }
                                    default:
                                        break;
                                }
                                using namespace Acts::detail::LineHelper;
                                const Acts::Intersection3D bsExtp = lineIntersect<3>(Amg::Vector3D::Zero(),
                                                                                     Amg::Vector3D::UnitZ(),
                                                                                     globPos, globDir);
                                const Amg::Vector3D closePoint = bsExtp.position();
                                const auto rawStIdx = Acts::toUnderlying(stIdx);
                                //coverity[NEGATIVE_RETURNS]
                                if (closePoint.perp() > m_seedMaxBsR.value().at(rawStIdx) * m_beamSpotR ||
                                    std::abs(closePoint.z()) >  m_seedMaxBsL.value().at(rawStIdx) * m_beamSpotL){
                                    ATH_MSG_DEBUG("fitSegmentSeed() - Reject parameters "<<Amg::toString(tangentSeedPos)
                                                  <<" + "<<Amg::toString(tangentSeedDir)
                                                  <<" as extrapolation to beamspot is too far "<<Amg::toString(closePoint)
                                                  <<", r: "<<closePoint.perp());
                                    return false;
                                }
                                return true;
                          }};
 
        const auto* seeder = patternSeed->msSector()->chamberIndex() == ChIndex::BEE ?
                             m_seederBEE.get() : m_seeder.get();
   
        if (m_visionTool.isEnabled()) { 
            PrimitiveVec seedLines{};
            State_t drawMe{seedState};
            while(auto s = seeder->nextSeed(cctx, drawMe)) {
                seedLines.push_back(drawLine(s->parameters, -Gaudi::Units::m, Gaudi::Units::m, kViolet));
            }
            seedLines.push_back(drawLabel(std::format("possible seeds: {:d}",  drawMe.nGenSeeds()), 0.2, 0.85, 14));
            m_visionTool->visualizeSeed(ctx, *patternSeed, std::format("pattern_{:}{:}{:}",
                patternSeed->msSector()->chamberIndex(),
                patternSeed->msSector()->side() ? 'A' : 'C' ,
                patternSeed->msSector()->sector()), std::move(seedLines));
        }
 
        ATH_MSG_VERBOSE("fitSegmentHits() - Start segment seed search");
        while (auto seed = seeder->nextSeed(cctx, seedState)) {
            ATH_MSG_VERBOSE("fitSegmentHits() - Found a seed. Try to fit the segment...");
            // Back convert the seed parameters to athena units
            seed->parameters[toUnderlying(ParamDefs::t0)] = ActsTrk::timeToAthena(seed->parameters[toUnderlying(ParamDefs::t0)]);
            
            auto segment = m_fitter->fitSegment(ctx, patternSeed, seed->parameters,
                                                locToGlob, std::move(seed->hits));
            if (segment) {
                segments.push_back(std::move(segment));
            }
        }
        ATH_MSG_VERBOSE("fitSegmentHits() - In total "<<segments.size()<<" segment were constructed ");
        return segments;
 
    }
   
    void SegmentFittingAlg::resolveAmbiguities(const ActsTrk::GeometryContext& gctx,
                                               SegmentVec_t& segmentCandidates) const {
        if (segmentCandidates.empty()) {
            return;
        }
        ATH_MSG_VERBOSE("resolveAmbiguities() - Resolve ambiguities amongst "
                        <<segmentCandidates.size()<<" segment candidates.");
        std::map<const MuonGMR4::SpectrometerSector*, SegmentVec_t,
                 MuonGMR4::MuonDetectorManager::MSEnvelopeSorter> candidatesPerChamber{};
        
        for (std::unique_ptr<Segment>& sortMe : segmentCandidates) {
            const MuonGMR4::SpectrometerSector* chamb = sortMe->msSector();
            candidatesPerChamber[chamb].push_back(std::move(sortMe));
        }
        segmentCandidates.clear();
        for (auto& [chamber, resolveMe] : candidatesPerChamber) {
            const std::size_t nBefore = resolveMe.size();
            SegmentVec_t resolvedSegments = m_ambiSolver->resolveAmbiguity(gctx, std::move(resolveMe));
            ATH_MSG_DEBUG("resolveAmbiguities()  - "<<resolvedSegments.size()<<"/"<<nBefore
                <<" segments survived ambiguity solving in "<<chamber->identString()<<".");
            segmentCandidates.insert(segmentCandidates.end(), 
                                     std::make_move_iterator(resolvedSegments.begin()),
                                     std::make_move_iterator(resolvedSegments.end()));
        }
        ATH_MSG_VERBOSE("Ambiguity solving done "<<segmentCandidates.size()<<" survived.");
    }
}