MuonFastSABuilderAlg.cxx

Go to the documentation of this file.

/*
  Copyright (C) 2002-2026 CERN for the benefit of the ATLAS collaboration
*/
 
#include "MuonFastSABuilderAlg.h"
 
#include "MuonTrackEvent/MsTrackSeed.h"
#include "MuonFastRecoEvent/FastRecoUtils.h"
#include "Acts/Utilities/VectorHelpers.hpp"
 
namespace {
    double inDeg(double angle) {
        return angle / Gaudi::Units::deg;
    }
    std::string print(const xAOD::Muon& muon) {
        std::ostringstream oss;
        using namespace MuonR4::FastReco;
        const double P {muon.pt() * std::cosh(muon.eta())};
        const double momUncertainty {std::sqrt(getQOverPCov(muon)*std::pow(P,4))};
        oss<<"Eta: "<<muon.eta()<<", Phi: "<<inDeg(muon.phi())
           <<", Pt: "<<muon.pt()/Gaudi::Units::GeV<<" GeV, Charge: "<<muon.charge()
           <<", P: "<<P/ Gaudi::Units::GeV<<" += "<<momUncertainty/Gaudi::Units::GeV
           <<" Gev"<<std::endl<<"Associated "<< *getParentPattern(muon);
        return oss.str();
    }
}
 
namespace MuonR4 {
using namespace Acts::UnitLiterals;
using namespace Muon::MuonStationIndex;
 
StatusCode MuonFastSABuilderAlg::initialize() {
    ATH_CHECK(m_outMuons.initialize());
    ATH_CHECK(m_qOverPCovKey.initialize());
    ATH_CHECK(m_inSegments.initialize());
    ATH_CHECK(m_seedingTool.retrieve());
    return StatusCode::SUCCESS;
}
 
StatusCode MuonFastSABuilderAlg::execute(const EventContext& ctx) const {
    ATH_MSG_VERBOSE(__func__<<"() Start building muon candidates");
 
    const xAOD::MuonSegmentContainer* inSegments{nullptr};
    ATH_CHECK(SG::get(inSegments, m_inSegments, ctx));
 
    MuonDataShip dataShip{};
    ATH_CHECK(dataShip.muonContainer.record(m_outMuons, ctx));
    dataShip.dec_qOverPCov.initialize(m_qOverPCovKey, ctx);

    auto firstSeg = inSegments->stdcont().begin();
    const auto lastSeg = inSegments->stdcont().end();
    while(firstSeg != lastSeg) {
        const GlobalPattern* pattern {FastReco::getParentPattern(**firstSeg)};
        assert(pattern);
 
        const auto nextToLastSeg {std::find_if(firstSeg, lastSeg,
            [pattern](const xAOD::MuonSegment* seg) {
                return FastReco::getParentPattern(*seg) != pattern;
            })};
        const std::span<const xAOD::MuonSegment* const> muonSegments {firstSeg, nextToLastSeg};
 
        xAOD::Muon* newMuon {buildMuonCandidate(ctx, *pattern, muonSegments, dataShip)};
 
        firstSeg = nextToLastSeg;
        if (!newMuon) {
            ATH_MSG_DEBUG(__func__<<"() No muon candidate could be built from pattern " << *pattern);
            continue;
        }
        std::vector<ElementLink<xAOD::MuonSegmentContainer>> segLinks{};
        for (const xAOD::MuonSegment* seg : muonSegments) {
            segLinks.emplace_back(*inSegments, seg->index());
        }
        newMuon->setMuonSegmentLinks(segLinks);
        ATH_MSG_DEBUG(__func__<<"() Built new muon candidate: " << print(*newMuon));
    }
    ATH_MSG_DEBUG("Written "<<dataShip.muonContainer->size()<<" FastMuonSA into StoreGate.");
    return StatusCode::SUCCESS;
}
 
xAOD::Muon* MuonFastSABuilderAlg::buildMuonCandidate(const EventContext& ctx,
                                                     const GlobalPattern& pattern,
                                                     std::span<const xAOD::MuonSegment* const> segments,
                                                     MuonDataShip& outMuonData) const {
 
    ATH_MSG_VERBOSE(__func__<<"() Starting muon building from "
        << segments.size()<<" segments for pattern " << std::endl << pattern);
 
    if (segments.size() < 2) {
        ATH_MSG_DEBUG(__func__<<"() Not enough muon segments to construct a candidate - abort.");
        return nullptr;
    }

    const std::vector<StIndex> stations {pattern.getStations()};
    const unsigned nBarrelStations = std::ranges::count_if(stations, 
        [](const StIndex station) { return isBarrel(station); });
    
    using enum MsTrackSeed::Location;
    MsTrackSeed trackSeed {nBarrelStations * 2u >= stations.size() ? Barrel : Endcap, 
                           pattern.expSector()};
    for (const xAOD::MuonSegment* seg : segments) {
        trackSeed.addSegment(seg);
    }
    ATH_MSG_VERBOSE(__func__<<"() Constructed track seed: " << trackSeed);
 
    Acts::Result<Acts::BoundTrackParameters> seedPars {
        m_seedingTool->estimateStartParameters(ctx, trackSeed)};
 
    if (!seedPars.ok()) {
        ATH_MSG_DEBUG(__func__<<"() Failed to estimate seed parameters for track seed: " << trackSeed);
        return nullptr;
    }
 
    xAOD::Muon* newMuon = outMuonData.muonContainer->push_back(std::make_unique<xAOD::Muon>());
    newMuon->setAuthor(xAOD::Muon::Author::MuidSA);
    newMuon->setMuonType(xAOD::Muon::MuonType::MuonStandAlone);
    newMuon->setCharge(seedPars->charge());
 
    newMuon->setP4(ActsTrk::energyToAthena(seedPars->transverseMomentum()), 
                   Acts::VectorHelpers::eta(*seedPars),seedPars->phi());
    
    std::optional<Acts::BoundMatrix> cov {seedPars->covariance()};
    float& qOverPCov {outMuonData.dec_qOverPCov(*newMuon)};
    if (!cov || (*cov)(Acts::eBoundQOverP, Acts::eBoundQOverP) < Acts::s_epsilon) {
        ATH_MSG_WARNING(__func__<<"() No covariance matrix available for seed parameters of track seed: " << trackSeed);
        qOverPCov = 0.;
        return newMuon;
    }
    static constexpr double conversionFactor {1 / Acts::square(ActsTrk::energyToAthena(1.))};
    qOverPCov = (*cov)(Acts::eBoundQOverP, Acts::eBoundQOverP) * conversionFactor;
    return newMuon;
}
}