PrepDataToSimHitAssocAlg.cxx

Go to the documentation of this file.

/*
  Copyright (C) 2002-2024 CERN for the benefit of the ATLAS collaboration
*/
#include "PrepDataToSimHitAssocAlg.h"
 
#include "StoreGate/ReadHandle.h"
#include "StoreGate/WriteDecorHandle.h"
#include "xAODMuonPrepData/UtilFunctions.h"
#include "xAODMuonPrepData/RpcMeasurement.h"
#include "xAODMuonPrepData/TgcStrip.h"
#include "xAODMuonPrepData/sTgcMeasurement.h"
#include "xAODMuonPrepData/MMCluster.h"
#include "xAODMuonViews/ChamberViewer.h"
namespace MuonR4{
    template <class ContainerType>
        StatusCode PrepDataToSimHitAssocAlg::retrieveContainer(const EventContext& ctx, 
                                                               const SG::ReadHandleKey<ContainerType>& key,
                                                               const ContainerType*& contToPush) const {
        contToPush = nullptr;
        if (key.empty()) {
            ATH_MSG_VERBOSE("No key has been parsed for object "<< typeid(ContainerType).name());
            return StatusCode::SUCCESS;
        }
        SG::ReadHandle readHandle{key, ctx};
        ATH_CHECK(readHandle.isPresent());
        contToPush = readHandle.cptr();
        return StatusCode::SUCCESS;
    }
 
    StatusCode PrepDataToSimHitAssocAlg::initialize() {
        ATH_CHECK(m_simHitsKey.initialize());
        ATH_CHECK(m_prdHitKey.initialize());
        ATH_CHECK(m_decorKey.initialize());
        ATH_CHECK(m_idHelperSvc.retrieve());
        ATH_CHECK(m_geoCtxKey.initialize());
        return StatusCode::SUCCESS;
    }    
    StatusCode PrepDataToSimHitAssocAlg::execute(const EventContext & ctx) const {
        const ActsGeometryContext* gctx{nullptr};
        const xAOD::MuonSimHitContainer* simHits{nullptr};
        const xAOD::UncalibratedMeasurementContainer* measurements{nullptr};
        ATH_CHECK(retrieveContainer(ctx, m_geoCtxKey, gctx));
        ATH_CHECK(retrieveContainer(ctx, m_simHitsKey, simHits));
        ATH_CHECK(retrieveContainer(ctx, m_prdHitKey, measurements));
        
        
        SG::WriteDecorHandle<xAOD::UncalibratedMeasurementContainer, LinkType> decorHandle{m_decorKey, ctx};
        for (const xAOD::UncalibratedMeasurement* measurement : *measurements){
            const xAOD::MuonSimHit* bestSimHit{nullptr};
 
            switch (measurement->type()) {
                case xAOD::UncalibMeasType::MdtDriftCircleType: {
                    const Identifier prdId{xAOD::identify(measurement)};
                    xAOD::MuonSimHitContainer::const_iterator mdt_matching = 
                        std::ranges::find_if(*simHits,[&prdId](const xAOD::MuonSimHit* hit){
                            return hit->identify() == prdId;
                        });
                    if (mdt_matching != simHits->end()){
                        bestSimHit =(*mdt_matching);
                    }
                    break;
                }
                case xAOD::UncalibMeasType::RpcStripType:
                case xAOD::UncalibMeasType::TgcStripType:
                case xAOD::UncalibMeasType::MMClusterType:
                case xAOD::UncalibMeasType::sTgcStripType: {
 
                    const Identifier prdId{xAOD::identify(measurement)};
                    const MuonGMR4::MuonReadoutElement* readOutEle = xAOD::readoutElement(measurement);
                    const Amg::Transform3D& locToGlob{readOutEle->localToGlobalTrans(*gctx, readOutEle->layerHash(prdId))};
                    
                    const Identifier gasGapId = m_idHelperSvc->gasGapId(prdId);
                    Amg::Vector3D locPos{Amg::Vector3D::Zero()};
                    if (measurement->numDimensions() == 1) {
                        locPos = measurement->localPosition<1>().x() * Amg::Vector3D::UnitX();
                    } else {
                        locPos.block<2,1>(0,0) = xAOD::toEigen(measurement->localPosition<2>());
                    }
                    double closestDistance{m_PullCutOff};
 
                    xAOD::ChamberViewer chambViewer{*simHits, m_idHelperSvc.get(), 
                                                    xAOD::ChamberView::Mode::DetElement};
 
                    if (!chambViewer.loadView(gasGapId)) {
                        break;
                    }
                    for ( const xAOD::MuonSimHit* simHit : chambViewer) {
                        if (gasGapId != m_idHelperSvc->gasGapId(simHit->identify())) {
                            continue; 
                        }
                        const IdentifierHash simLayHash{readOutEle->layerHash(simHit->identify())};
                        const Amg::Transform3D globToLoc{readOutEle->globalToLocalTrans(*gctx, simLayHash) *locToGlob};
                        const Amg::Vector3D prdPos = globToLoc * locPos;
                        double dist{0.};
                        if (measurement->numDimensions() == 1) {
                            dist = std::abs(prdPos.x() - simHit->localPosition().x()) 
                                 / std::sqrt(measurement->localCovariance<1>()(0,0));
                        } else{
                            const Amg::Vector2D diff = (prdPos - xAOD::toEigen(simHit->localPosition())).block<2,1>(0,0);
                            dist = std::sqrt(diff.dot(xAOD::toEigen(measurement->localCovariance<2>()).inverse() * diff));
                        }
                        if (dist < closestDistance) {
                            closestDistance = dist;
                            bestSimHit = simHit;
                        }
                    }                    
                    break;
                } default: {
                    ATH_MSG_FATAL("Non muon measurement is parsed");
                    return StatusCode::FAILURE;
                }
            }
            if (!bestSimHit) {
                continue;
            }
            decorHandle(*measurement) = LinkType{*simHits, bestSimHit->index()};
        }
        return StatusCode::SUCCESS;
    }
}