xAODSegmentCnvTool.cxx

Go to the documentation of this file.

/*
  Copyright (C) 2002-2026 CERN for the benefit of the ATLAS collaboration
*/
 
#include "xAODSegmentCnvTool.h"
 
#include "EventPrimitives/EventPrimitivesHelpers.h"
#include "xAODMuonPrepData/UtilFunctions.h"
#include "MuonTrackEvent/TrackingHelpers.h"
#include "MuonPatternEvent/MuonHoughDefs.h"
 
#include "Acts/Surfaces/StrawSurface.hpp"
#include "ActsGeoUtils/SurfacePlacement.h"
#include "Acts/Geometry/TrackingGeometry.hpp"
#include "Acts/Surfaces/PlaneSurface.hpp"
#include "Acts/Definitions/Units.hpp"
 
namespace MuonR4{
 
   using namespace Acts::UnitLiterals;
 
 
    StatusCode xAODSegmentCnvTool::initialize() {
        ATH_CHECK(m_idHelperSvc.retrieve());
        ATH_CHECK(m_ctxProvider.initialize());
        if (m_estimateHoles) {
            ATH_CHECK(m_trackingGeometrySvc.retrieve());
        }
        ATH_CHECK(m_extrapolationTool.retrieve(EnableTool{m_estimateHoles}));
        return StatusCode::SUCCESS;
    }
 
    xAOD::MuonSegment* xAODSegmentCnvTool::convertSegment(const EventContext& ctx,
                                                      const MuonR4::Segment& inSegment,
                                                      DataShip& ship) const {
 
        const Acts::GeometryContext tgContext = m_ctxProvider.getGeometryContext(ctx);
        const MuonGMR4::SpectrometerSector* sector = inSegment.msSector();
 
        xAOD::MuonSegment* convertedSeg = ship.segmentContainer->push_back(std::make_unique<xAOD::MuonSegment>());
        
        const Amg::Vector3D& pos{inSegment.position()};
        const Amg::Vector3D& dir{inSegment.direction()};
        convertedSeg->setPosition(pos.x(), pos.y(), pos.z());
        convertedSeg->setDirection(dir.x(), dir.y(), dir.z());
 
        convertedSeg->setIdentifier(sector->sector(), sector->chamberIndex(),
                                    sector->side(), inSegment.technology());
        convertedSeg->setFitQuality(inSegment.chi2(), inSegment.nDoF());
 
        using enum SegmentFit::ParamDefs;
        convertedSeg->setT0Error(inSegment.segementT0(), 
                                 Amg::error(inSegment.covariance(), Acts::toUnderlying(t0)));
        
        if (ship.dec_localSegPars) {
            using SegPars_t = xAOD::PosAccessor<Acts::toUnderlying(SegmentFit::ParamDefs::nPars)>::element_type;
            SegPars_t& localPars{ship.dec_localSegPars(*convertedSeg)};
        
            const Amg::Transform3D globToLoc{sector->globalToLocalTransform(tgContext)};            
            const Amg::Vector3D locPos{globToLoc * pos};
            const Amg::Vector3D locDir{globToLoc.linear() * dir};
 
            localPars[Acts::toUnderlying(x0)] = locPos.x();
            localPars[Acts::toUnderlying(y0)] = locPos.y();
            localPars[Acts::toUnderlying(theta)] = locDir.theta();
            localPars[Acts::toUnderlying(phi)] = locDir.phi();
            localPars[Acts::toUnderlying(t0)] = inSegment.segementT0();
        }
 
        if (ship.dec_localSegCov) {
            using SegCov_t = xAOD::PosAccessor<Acts::sumUpToN(Acts::toUnderlying(SegmentFit::ParamDefs::nPars))>::element_type;
            SegCov_t& localCov{ship.dec_localSegCov(*convertedSeg)};
            constexpr std::size_t n = Acts::toUnderlying(SegmentFit::ParamDefs::nPars);
            for (std::size_t p = 1; p < n; ++p) {
                for (std::size_t p1 = 0 ; p1 <=p;++p1) {
                    localCov[Acts::vecIdxFromSymMat<n>(p,p1)] = inSegment.covariance()(p, p1);
                }
            }
        }
 
        if (linkMeasurements(tgContext, inSegment, *convertedSeg, ship) != StatusCode::SUCCESS) {
            ATH_MSG_ERROR("Failed to link measurements");
            ship.segmentContainer->pop_back();
            return nullptr;
        }
        evaluateSummary(ctx, *convertedSeg);
        return convertedSeg;
    }
 
    StatusCode xAODSegmentCnvTool::linkMeasurements(const Acts::GeometryContext& tgContext,
                                                const Segment& inSegment,
                                                xAOD::MuonSegment& targetSegment,
                                                DataShip& ship) const {
        if (!ship.dec_prdLinks || !ship.dec_prdStates || !ship.prdCombContainer.hasHandle()) {
            ATH_MSG_DEBUG("Decorators for the prd links and states are not initialized");
            return StatusCode::SUCCESS;
        }
 
        using State = CalibratedSpacePoint::State;
        // Cache all measurements that can be combined to two measurements in a single gas gap
        std::vector< std::tuple<const xAOD::MuonMeasurement*, State, std::size_t>> combineMap{};
        std::vector< std::tuple<const xAOD::UncalibratedMeasurement*, State, std::size_t>> linkMap{};
 
        using PrdLinkVec_t = std::vector<ElementLink<xAOD::UncalibratedMeasurementContainer>>;
        PrdLinkVec_t& links = ship.dec_prdLinks(targetSegment);
        std::vector<char>& linkStates = ship.dec_prdStates(targetSegment);
        links.reserve(2*inSegment.measurements().size());
        linkStates.reserve(2*inSegment.measurements().size());

        auto combine = [&](const xAOD::MuonMeasurement* m1, 
                          const xAOD::MuonMeasurement* m2) {
            auto cmbMeas = ship.prdCombContainer->push_back(std::make_unique<xAOD::CombinedMuonStrip>());
 
            cmbMeas->setPrimaryStrip(m1);
            cmbMeas->setSecondaryStrip(m2);
            const auto [locPos, locCov] = xAOD::positionAndCovariance(m1, m2);
            cmbMeas->localCovariance<2>() = xAOD::toStorage(locCov);
            cmbMeas->localPosition<2>() = xAOD::toStorage(locPos);
            const Identifier id1{m1->identify()}, id2{m2->identify()};
            ATH_MSG_VERBOSE(__func__<<"() "<<__LINE__<<" - Combine "<<m_idHelperSvc->toString(id1) 
                            <<" & "<<m_idHelperSvc->toString(id2));
            if ((m1->type() != xAOD::UncalibMeasType::sTgcStripType || 
                m_idHelperSvc->stgcIdHelper().channelType(id1) == 
                m_idHelperSvc->stgcIdHelper().channelType(id2))&&
                m_idHelperSvc->measuresPhi(id1) == m_idHelperSvc->measuresPhi(id2)) {
                THROW_EXCEPTION("Cannot combine "<<m_idHelperSvc->toString(id1)
                                <<" & "<<m_idHelperSvc->toString(id2));
            }
            return cmbMeas;
        };
        // Loop over the measurements
        for (const auto [segIdx, meas] : Acts::enumerate(inSegment.measurements())) {
            const SpacePoint* sp = meas->spacePoint();
            if (!sp) {
                if (!ship.beamSpotMeasCreator.ok()) {
                    continue;
                }
                // Up to now, there's no variety on the beamspot across the segments
                if (!ship.beamSpot) {
                    if (!ship.beamSpotMeasCreator.ok()) {
                        ATH_MSG_ERROR("Cannot create a beamspot measurement");
                        return StatusCode::FAILURE;
                    }
 
                    const Amg::Vector3D beamSpot = inSegment.msSector()->localToGlobalTransform(tgContext) *
                                                    meas->localPosition();
                    AmgSymMatrix(2) covariance{AmgSymMatrix(2)::Identity()};
                    using CovIdx = SpacePoint::CovIdx;
                    using ProjectorType = xAOD::AuxiliaryMeasurement::ProjectorType;
                    covariance(0,0) = meas->covariance()[Acts::toUnderlying(CovIdx::etaCov)];
                    covariance(1,1) = meas->covariance()[Acts::toUnderlying(CovIdx::phiCov)];
                    auto surf = Acts::Surface::makeShared<Acts::StrawSurface>(Amg::getTranslate3D(beamSpot),
                                                std::make_shared<Acts::LineBounds>(std::sqrt(covariance(0,0)), 20._m));
 
                    ship.beamSpot = ship.beamSpotMeasCreator->newMeasurement<2>(surf, ProjectorType::e2DimNoTime, covariance);
                    ATH_MSG_DEBUG(__func__<<"() "<<__LINE__<<" - Created beamspot measurement "<<(*meas)<<", "
                                    <<surf->toString(tgContext));
                }
                linkMap.emplace_back(ship.beamSpot, meas->fitState(), segIdx);
                continue;
            }
            switch (sp->type()) {
                using enum xAOD::UncalibMeasType;
                // Mdt &  micromegas are never combined
                case MdtDriftCircleType:
                case MMClusterType: {
                    linkMap.emplace_back(sp->primaryMeasurement(), meas->fitState(), segIdx);                        
                    break;
                } case RpcStripType:
                  case TgcStripType:
                  case sTgcStripType: {
                    if (sp->primaryMeasurement() && sp->secondaryMeasurement()) {
                        if (sp->primaryMeasurement() != sp->secondaryMeasurement()) {
                            linkMap.emplace_back(combine(sp->primaryMeasurement(), 
                                                         sp->secondaryMeasurement()),
                                                         meas->fitState(), segIdx);
                        } else {  // BI - RPC measurements
                            linkMap.emplace_back(sp->primaryMeasurement(), meas->fitState(), segIdx);
                        }
                    } else {
                        ATH_MSG_VERBOSE(__func__<<"() "<<__LINE__<<" - Append for later combination "<<(*meas));
                        combineMap.emplace_back(sp->primaryMeasurement(), meas->fitState(), segIdx);
                    }
                    break;
                } default:
                    break;
            }
        }
    
        // Finally we need to check whether there're measurements left to combine
        for (std::size_t cmbIdx = 0; cmbIdx < combineMap.size(); ++cmbIdx){
            const xAOD::MuonMeasurement* m1{std::get<0>(combineMap[cmbIdx])};
            const State s1{std::get<1>(combineMap[cmbIdx])};
            const std::size_t segIdx1{std::get<2>(combineMap[cmbIdx])};
            ATH_MSG_VERBOSE(__func__<<"() "<<__LINE__<<" - Find another measurement to combine with "
                            <<m_idHelperSvc->toString(m1->identify()));
            if (cmbIdx +1 < combineMap.size()) {
                const xAOD::MuonMeasurement* m2{std::get<0>(combineMap[cmbIdx +1])};
                const State s2{std::get<1>(combineMap[cmbIdx+1])};
                ATH_MSG_VERBOSE(__func__<<"() "<<__LINE__<<" - Check whether "<<m_idHelperSvc->toString(m2->identify())
                                <<" is a good candidate");
                // Ensure that they point to the same detector element, are within the same
                // layer and have also the same state
                if (m1->type() == m2->type() &&  m1->identifierHash() == m2->identifierHash() && 
                    m1->layerHash() == m2->layerHash() && s1 == s2) {
                    ATH_MSG_VERBOSE(__func__<<"() "<<__LINE__<<" - They match");
                    if (m1->measuresPhi()) {
                        linkMap.emplace_back(combine(m2, m1), s2, segIdx1);
                    } else {
                        linkMap.emplace_back(combine(m1, m2), s1, segIdx1);
                    }
                    ++cmbIdx; // skip the next measurement as it's absorbed here
                    continue;
                }
            }
            ATH_MSG_VERBOSE("No match found");
            linkMap.emplace_back(m1, s1, segIdx1);
        }
 
        std::ranges::sort(linkMap, [](const auto& a, const auto& b){
            return std::get<2>(a) < std::get<2>(b);
        });
 
        for (const auto& [prd, state, segIdx]: linkMap) {
            links.emplace_back(*static_cast<const xAOD::UncalibratedMeasurementContainer*>(prd->container()), 
                                prd->index());
            linkStates.emplace_back(Acts::toUnderlying(state));
        }
        return StatusCode::SUCCESS;
    }
 
    void xAODSegmentCnvTool::evaluateSummary(const EventContext& ctx, 
                                         xAOD::MuonSegment& segment) const {

        Counter hits{}, outliers{}, holes{};
        std::unordered_set<Identifier> crossedSurfaces{};
        const std::size_t nMeas = nMeasurements(segment);
 
        const Acts::Surface* startSurface{}, *lastSurface{};
        const Acts::TrackingGeometry* trackingGeo = m_estimateHoles ? m_trackingGeometrySvc->trackingGeometry().get() : nullptr;
 
        ATH_MSG_DEBUG(__func__<<"() "<<__LINE__<<" - Loop over "<<nMeas<<" measurements.");
        for (std::size_t m = 0 ; m < nMeas; ++m) {
            const xAOD::UncalibratedMeasurement* meas = getMeasurement(segment, m);
            // Skip auxiliary constraints from summary, such as the beam spot (technically this should not happen because the beam spot is not considered as measurement here, but just in case )
            if (meas->type() == xAOD::UncalibMeasType::Other) {
                ATH_MSG_VERBOSE(__func__ << "() " << __LINE__
                                << " - Skip auxiliary measurement in segment summary");
                continue;
            }
 
            const bool isOutlier = isOutlierMeasurement(segment, m);
            Counter& increment = {isOutlier ? outliers: hits};
 
            increment.precision  += xAOD::isPrecisionHit(meas);
            const auto* muonMeas = static_cast<const xAOD::MuonMeasurement*>(meas);
 
            increment.triggerPhi += (muonMeas->measuresPhi()  ||
                                    //RPC BI / MDT twin or sTGC strip + pad hits
                                    muonMeas->numDimensions() == 2);
            increment.triggerEta += !xAOD::isPrecisionHit(muonMeas) && ( 
                                        // 0D for combined strip + X measurements
                                        (muonMeas->numDimensions() == 0 ) || 
                                        // 1D strip/wire/pad (RPC,TGC, sTGC) cases
                                        (muonMeas->numDimensions() == 1 && !muonMeas->measuresPhi() ) || 
                                        //RPC BI for which measuresPhi = 0 as the phi is extracted from time information
                                        (muonMeas->numDimensions() == 2));
            // Count measurement holes of the trigger hits
            if (!isOutlier && meas->numDimensions() == 1) {
                if (meas->type() == xAOD::UncalibMeasType::RpcStripType) {
                  //coverity[FORWARD_NULL]
                    const auto* re = static_cast<const MuonGMR4::RpcReadoutElement*>(muonMeas->readoutElement());
                    if (!muonMeas->measuresPhi() && re->nPhiStrips()) {
                        ++holes.triggerPhi;
                    } else if (muonMeas->measuresPhi()) {
                        ++holes.triggerEta;
                    }
                } else if (meas->type() == xAOD::UncalibMeasType::TgcStripType) {
                    //coverity[FORWARD_NULL]
                    const auto* re = static_cast<const MuonGMR4::TgcReadoutElement*>(muonMeas->readoutElement());
                        if (!muonMeas->measuresPhi() && re->numStrips(muonMeas->layerHash())) {
                        ++holes.triggerPhi;
                    } else if (muonMeas->measuresPhi() && re->numWireGangs(muonMeas->layerHash())) {
                        ++holes.triggerEta;
                    }
                } else if (meas->type() == xAOD::UncalibMeasType::sTgcStripType)  {
                    if (!xAOD::isPrecisionHit(meas)) {
                        ++holes.precision;
                    } else {
                        ++holes.triggerPhi;
                    }
                }
            }
            
            if (!m_estimateHoles) {
                continue;
            }   
            const Acts::Surface& surface = xAOD::muonSurface(meas);
            if (!muonMeas){
                continue;
            }
            crossedSurfaces.insert( muonMeas->type() != xAOD::UncalibMeasType::MdtDriftCircleType ? 
                                    m_idHelperSvc->gasGapId(muonMeas->identify()) :
                                    muonMeas->identify());
            const auto* volume = MuonGMR4::highestAlignable(trackingGeo->findVolume(volumeId(surface)));
            assert(volume != nullptr);
            if (!startSurface) {
                startSurface = MuonGMR4::bottomBoundary(*volume);
            } else if (m +1 == nMeas) {
                lastSurface = MuonGMR4::topBoundary(*volume);
            }
        }
        if (m_estimateHoles) {
            const Acts::GeometryContext tgContext = m_ctxProvider.getGeometryContext(ctx);
            std::shared_ptr<const Acts::Surface> startSurfPtr = startSurface->getSharedPtr();
            auto atSurface = startSurface->intersect(tgContext, segment.position(), segment.direction(),
                                                      Acts::BoundaryTolerance::Infinite()).closest();
 
            const auto locPos = startSurface->globalToLocal(tgContext, atSurface.position(), segment.direction());
            if (!locPos.ok()) {
                ATH_MSG_WARNING(__func__<<"() "<<__LINE__<<" - Intersection "<<Amg::toString(atSurface.position())
                    <<" is not on surface "<<startSurface->geometryId());
                return;
            }
            if (!startSurface->insideBounds(*locPos)) {
                const auto* volume = MuonGMR4::highestAlignable(trackingGeo->highestTrackingVolume()->
                                                      lowestTrackingVolume(tgContext, atSurface.position()));
                ATH_MSG_DEBUG(__func__<<"() "<<__LINE__<<" - Intersection "<<Amg::toString(atSurface.position())
                    <<" is outside surface: "<<startSurface->geometryId()<<", "<<startSurface->bounds()<<". "
                    <<" Switch to volume "<<volume->geometryId()<<", bounds: "<<volume->volumeBounds());
 
                if (volume->isAlignable()) {
                    startSurfPtr = MuonGMR4::bottomBoundary(*volume)->getSharedPtr();
                    atSurface = startSurfPtr->intersect(tgContext, segment.position(), segment.direction(),
                                                        Acts::BoundaryTolerance::Infinite()).closest();
                } else {
                    startSurfPtr = Acts::Surface::makeShared<Acts::PlaneSurface>(startSurfPtr->localToGlobalTransform(tgContext));
                }
            }
            auto startPars = Acts::BoundTrackParameters::create(tgContext, startSurfPtr,
                                                                ActsTrk::convertPosToActs(atSurface.position()),
                                                                segment.direction(), 1./ 5._TeV, std::nullopt,
                                                                Acts::ParticleHypothesis::muon());
            if (startPars.ok()) {
                findHoles(ctx, *startPars, lastSurface, crossedSurfaces, holes);
            } else {
                ATH_MSG_WARNING(__func__<<"() "<<__LINE__<<" - Start parameters not defined.");
            }
        }
        segment.setNHits(hits.precision, hits.triggerPhi, hits.triggerEta);
        segment.setNOutliers(outliers.precision, outliers.triggerPhi, outliers.triggerEta);
        segment.setNHoles(holes.precision, holes.triggerPhi, holes.triggerEta);
        ATH_MSG_VERBOSE(__func__<<"() "<<__LINE__<<" Updated summary: "<<hits);
    }
    void xAODSegmentCnvTool::findHoles(const EventContext& ctx,
                                      const Acts::BoundTrackParameters& startPars,
                                      const Acts::Surface* target,
                                      const std::unordered_set<Identifier>& layersWithHits,
                                      Counter& holeCounter) const {
        
 
        using SurfaceRecordOptions = ActsTrk::IExtrapolationTool::SurfaceRecordOptions;
        SurfaceRecordOptions propOpts{target, m_extraHolePath};
        propOpts.recordMaterial = false;
        propOpts.recordPassive = false;
        propOpts.recordSensitive = true;
 
        
        auto propResult = m_extrapolationTool->propagateAndRecord(ctx, startPars, propOpts);
        if (!propResult.ok()) {
            return;
        }
 
        /* Loop over the track record to filter out the holes*/
        for (Acts::BoundTrackParameters& record : (*propResult)) {
            const auto* placement = dynamic_cast<const ActsTrk::SurfacePlacement*>(record.referenceSurface().surfacePlacement());
            assert(placement != nullptr);
            const auto* detEl = static_cast<const MuonGMR4::MuonReadoutElement*>(placement->detectorElement());
            using enum ActsTrk::DetectorType;
            if (layersWithHits.count(placement->identify())) {
                continue;
            } 
            if(std::ranges::none_of(layersWithHits,[&](const Identifier& recorded){
                return m_idHelperSvc->chamberIndex(recorded) == detEl->chamberIndex() &&
                        m_idHelperSvc->stationPhi(recorded) == detEl->stationPhi();
            })) {
                ATH_MSG_DEBUG(__func__<<"() "<<__LINE__<<" - Surface "
                    <<m_idHelperSvc->toStringGasGap(placement->identify())<<" is not in the same sector");
                continue;
            }
            switch (placement->detectorType()) {
                case Mm:
                case Mdt: {
                    ATH_MSG_VERBOSE(__func__<<"() "<<__LINE__<<" - Precision hole detected "
                                   <<m_idHelperSvc->toString(placement->identify()));
                    ++holeCounter.precision;
                    break;
                } case Rpc: {
                    const auto* re = static_cast<const MuonGMR4::RpcReadoutElement*>(placement->detectorElement());
                    ATH_MSG_VERBOSE(__func__<<"() "<<__LINE__<<" - Rpc hole detected "
                                <<m_idHelperSvc->toStringGasGap(placement->identify())<<".");
                    ++holeCounter.triggerEta;
                    holeCounter.triggerPhi += (re->nPhiStrips() > 0);
                    break;
                } case Tgc: {
                    const auto* re = static_cast<const MuonGMR4::TgcReadoutElement*>(placement->detectorElement());
                    holeCounter.triggerEta += (re->numWireGangs(placement->hash())>0);
                    holeCounter.triggerPhi += (re->numStrips(placement->hash())>0);
                    ATH_MSG_VERBOSE(__func__<<"() "<<__LINE__<<" - Tgc hole detected "
                                <<m_idHelperSvc->toStringGasGap(placement->identify())<<".");
                    break;
                } case sTgc: {
                    ATH_MSG_VERBOSE(__func__<<"() "<<__LINE__<<" - Precision hole detected "
                                   <<m_idHelperSvc->toString(placement->identify()));
                    ++holeCounter.precision;
                    ++holeCounter.triggerPhi;
                    break;
                } default:
                    ATH_MSG_WARNING(__func__<<"() "<<__LINE__<<" - Unknown detector type "
                                    <<placement->detectorType());
                    break;
            }
        }
    }
}