MuonLayerHoughTool.h

Go to the documentation of this file.

/*
  Copyright (C) 2002-2025 CERN for the benefit of the ATLAS collaboration
*/
 
#ifndef MUONHOUGHPATTERNTOOLS_MUONLAYERHOUGHTOOL_H
#define MUONHOUGHPATTERNTOOLS_MUONLAYERHOUGHTOOL_H
 
#include <set>
 
#include "AthenaBaseComps/AthAlgTool.h"
#include "FourMomUtils/xAODP4Helpers.h"
#include "GaudiKernel/ServiceHandle.h"
#include "GaudiKernel/ToolHandle.h"
#include "GeoPrimitives/GeoPrimitives.h"
#include "MuonClusterization/TgcHitClustering.h"
#include "MuonDetDescrUtils/MuonSectorMapping.h"
#include "MuonIdHelpers/IMuonIdHelperSvc.h"
#include "MuonLayerHough/MuonLayerHough.h"
#include "MuonLayerHough/MuonLayerHoughSelector.h"
#include "MuonLayerHough/MuonPhiLayerHough.h"
#include "MuonLayerHough/MuonRegionHough.h"
#include "MuonPattern/MuonPatternCombinationCollection.h"
#include "MuonPrepRawData/MMPrepDataContainer.h"
#include "MuonPrepRawData/MuonPrepDataContainer.h"
#include "MuonPrepRawData/sTgcPrepDataContainer.h"
#include "MuonRecHelperTools/MuonEDMPrinterTool.h"
#include "MuonRecToolInterfaces/HoughDataPerSec.h"
#include "MuonRecToolInterfaces/IMuonHoughPatternFinderTool.h"
#include "TFile.h"
#include "TTree.h"
#include "TrkTruthData/PRD_MultiTruthCollection.h"
#include "xAODMuon/MuonSegmentContainer.h"
#include "xAODTruth/TruthParticleContainer.h"
#include "MuonReadoutGeometry/MuonDetectorManager.h"
namespace Trk {
    class PrepRawData;
}
 
 
namespace MuonHough {
    class HitDebugInfo;
}
 
namespace Muon {
 
    class MuonLayerHoughTool :  public extends<AthAlgTool,  IMuonHoughPatternFinderTool> {
    public:
        using HashVec = std::vector<IdentifierHash>;
        using RegionHashVec = std::vector<HashVec>;
        using TechnologyRegionHashVec = std::vector<RegionHashVec>;
        using DetRegIdx = MuonStationIndex::DetectorRegionIndex;
 
        struct CollectionsPerSector {
            int sector{-1};
            TechnologyRegionHashVec technologyRegionHashVecs;
        };
        typedef std::vector<CollectionsPerSector> CollectionsPerSectorVec;
 
        using HitVec = HoughDataPerSec::HitVec;
        using RegionHitVec = HoughDataPerSec::RegionHitVec;
        using PhiHitVec = HoughDataPerSec::PhiHitVec;
        using RegionPhiHitVec = HoughDataPerSec::RegionPhiHitVec;
        using MaximumVec = HoughDataPerSec::MaximumVec;
        using PhiMaximumVec = HoughDataPerSec::PhiMaximumVec;
        using MaximumAssociationMap = HoughDataPerSec::MaximumAssociationMap;
        using RegionMaximumVec = HoughDataPerSec::RegionMaximumVec;
        using RegionPhiMaximumVec = HoughDataPerSec::RegionPhiMaximumVec;
        using HoughDataPerSector = HoughDataPerSec;
 
        class Road {
        public:
            Road(std::shared_ptr<MuonHough::MuonLayerHough::Maximum> seed_) : 
                seed(seed_) { add(seed_); }
            Road() = default;
            using DetRegIdx = MuonStationIndex::DetectorRegionIndex;
            DetRegIdx neighbouringRegion{DetRegIdx::DetectorRegionUnknown};
            int neighbouringSector{-1};
            std::shared_ptr<MuonHough::MuonLayerHough::Maximum> seed{nullptr};
            void add(std::shared_ptr<MuonHough::MuonLayerHough::Maximum> max) {
                maxima.emplace_back(max);
                maximumSet.insert(max);
            }
            void add(std::shared_ptr<MuonHough::MuonPhiLayerHough::Maximum> max) { phiMaxima.emplace_back(max); }
            MaximumVec maxima{};
            PhiMaximumVec phiMaxima{};
            
            std::set<std::shared_ptr<MuonHough::MuonLayerHough::Maximum>> maximumSet{};
 
            std::vector<MuonHough::MuonPhiLayerHough::Maximum> mergedPhiMaxima{};
        };
 
        using base_class::base_class;
        virtual ~MuonLayerHoughTool() = default;
 
        virtual StatusCode initialize() override;
 
        virtual std::pair<std::unique_ptr<MuonPatternCombinationCollection>, std::unique_ptr<HoughDataPerSectorVec>> find(
            const MdtPrepDataContainer* mdtCont, const CscPrepDataContainer* cscCols, const TgcPrepDataContainer* tgcCont,
            const RpcPrepDataContainer* rpcCont, const sTgcPrepDataContainer* stgcCont, const MMPrepDataContainer* mmCont,
            const EventContext& ctx) const override;
 
        virtual std::pair<std::unique_ptr<MuonPatternCombinationCollection>, std::unique_ptr<HoughDataPerSectorVec>> find(
            const std::vector<const MdtPrepDataCollection*>& mdtCols, const std::vector<const CscPrepDataCollection*>& cscCols,
            const std::vector<const TgcPrepDataCollection*>& tgcCols, const std::vector<const RpcPrepDataCollection*>& rpcCols,
            const MuonSegmentCombinationCollection*, const EventContext& ctx) const override;
 
    private:
        using TgcEdge = TgcClusterObj3D::Edge;
 
        void getSectors(const Amg::Vector3D& pos, std::vector<int>& sectors) const;
        void getSectors(const TgcClusterObj3D& tgc, std::vector<int>& sectors) const;
 
        double rCor(const Amg::Vector3D& pos, const Identifier& id) const;
        double rCor(const MuonCluster& rpc) const;
        double rCor(const MdtPrepData& mdt) const;
        double rCor(const TgcClusterObj3D& tgc, const TgcEdge val, int sector) const;
 
        int sublay(const Identifier& id, float z = 0) const;  // the z value is only used for the tgcs
 
        struct State {
            MaximumVec seedMaxima{};  // Does not own the contained objects, they're just references to objects stored in houghDataPerSectorVec.
            std::unique_ptr<HoughDataPerSectorVec> houghDataPerSectorVec{std::make_unique<HoughDataPerSectorVec>()};
            std::set<Identifier> truthHits;
            std::set<Identifier> foundTruthHits;
            std::set<Identifier> outputTruthHits;
        };
 
        std::pair<std::unique_ptr<MuonPatternCombinationCollection>, std::unique_ptr<HoughDataPerSectorVec>> analyse(State& state) const;
 
        void fillHitsPerSector(const EventContext& ctx, State& state, const int sector, const CollectionsPerSector& hashes,
                               const MdtPrepDataContainer* mdtCont, const CscPrepDataContainer* cscCont,
                               const TgcPrepDataContainer* tgcCont, const RpcPrepDataContainer* rpcCont,
                               const sTgcPrepDataContainer* stgcCont, const MMPrepDataContainer* mmCont) const;
 
        void fill(const EventContext& ctx, std::set<Identifier>& truthHits, const MdtPrepDataCollection& mdts, HitVec& hits) const;
        void fill(const EventContext& ctx, std::set<Identifier>& truthHits,
                  std::vector<std::unique_ptr<TgcHitClusteringObj>>& tgcClusteringObjs, const TgcPrepDataCollection& tgcs, HitVec& hits,
                  PhiHitVec& phiHits, int sector) const;
        void fill(const EventContext& ctx, std::set<Identifier>& truthHits, const RpcPrepDataCollection& rpcs, HitVec& hits,
                  PhiHitVec& phiHits) const;
        void fill(const EventContext& ctx, std::set<Identifier>& truthHits, const MMPrepDataCollection& mdts, HitVec& hits) const;
        void fill(const EventContext& ctx, std::set<Identifier>& truthHits, const sTgcPrepDataCollection& stgcs, HitVec& hits,
                  PhiHitVec& phiHits, int sector) const;
 
        void fill(const EventContext& ctx, std::set<Identifier>& truthHits, const CscPrepDataCollection& cscs, HitVec& hits,
                  PhiHitVec& phiHits) const;
 
        bool findMaxima(MaximumVec& seedMaxima,
                        MuonHough::MuonLayerHough& hough, HitVec& hits, MaximumVec& maxima) const;
        bool findMaxima(MuonHough::MuonPhiLayerHough& hough,
                        PhiHitVec& hits, PhiMaximumVec& maxima, int sector) const;
 
        void associateMaximaToPhiMaxima(MuonStationIndex::DetectorRegionIndex region, HoughDataPerSector& houghData,
                                        std::map<MuonHough::MuonPhiLayerHough::Maximum*, MaximumVec>& phiEtaAssociations,
                                        std::vector<MaximumVec>& unassEtaMaxima) const;
 
        void associateMaximaInNeighbouringSectors(HoughDataPerSector& houghData,
                                                  std::vector<HoughDataPerSector>& houghDataPerSectorVec) const;
 
        void extendSeed(MuonHough::MuonDetectorHough& detectorHoughTransforms, 
                        Road& road, HoughDataPerSector& sectorData) const;  // const;
        void associatePhiMaxima(Road& road, PhiMaximumVec& phiMaxima) const;
 
        double combinedPeakheight(double ph, double ph1, double ph2, double phn, double rot, int layer, int /*region*/) const;
        void createPatternCombinations(std::vector<MaximumVec>& maxima, MuonPatternCombinationCollection& patternCombis) const;
 
        void createPatternCombinations(std::map<MuonHough::MuonPhiLayerHough::Maximum*, MaximumVec>& phiEtaAssociations,
                                       MuonPatternCombinationCollection& patternCombis) const;
 
        void insertHash(const IdentifierHash& hash, const Identifier& id) const;
        void insertHash(int sector, const IdentifierHash& hash, const Identifier& id) const;
 
        void matchTruth(std::set<Identifier>& truthHits, const PRD_MultiTruthCollection& truthCol, const Identifier& id,
                        MuonHough::HitDebugInfo& debug) const;
        
        void initializeSectorMapping(const EventContext& ctx) const;
        void printTruthSummary(std::set<Identifier>& truth, std::set<Identifier>& found) const;
 
        void buildRoads(MaximumVec& seedMaxima, MuonHough::MuonDetectorHough& detectorHoughTransforms, 
                        std::unique_ptr<HoughDataPerSectorVec>& houghDataPerSectorVec,
                        std::vector<Road>& roads) const;
        void mergePhiMaxima(Road& road) const;
 
        Gaudi::Property<bool> m_useSeeds{this, "UseSeeds", true};
        
        SG::ReadCondHandleKey<MuonGM::MuonDetectorManager> m_muonManagerKey{this, "MuonManagerKey", "MuonDetectorManager", "MuonManager ReadKey for IOV Range intersection"};
 
 
        ServiceHandle<IMuonIdHelperSvc> m_idHelperSvc{this, "MuonIdHelperSvc", "Muon::MuonIdHelperSvc/MuonIdHelperSvc"};
        PublicToolHandle<MuonEDMPrinterTool> m_printer{this, "printerTool", "Muon::MuonEDMPrinterTool/MuonEDMPrinterTool"};
 
        std::vector<MuonHough::MuonLayerHoughSelector> m_selectors;
        std::vector<MuonHough::MuonLayerHoughSelector> m_selectorsLoose;
 
        SG::ReadHandleKeyArray<PRD_MultiTruthCollection> m_truthNames{this, "TruthNames", {}};
 
        Gaudi::Property<bool> m_useRpcTimeVeto{this, "RpcTimeVeto", false};
        Gaudi::Property<bool> m_requireTriggerConfirmationNSW{this, "TriggerConfirmationNSW", false};
        Gaudi::Property<bool> m_onlyUseCurrentBunch{this, "OnlyUseCurrentBunch", false};
        Gaudi::Property<bool> m_debugHough{this, "DebugHough", false};
        Gaudi::Property<bool> m_doParabolicExtrapolation{this, "DoParabolicExtrapolation",
                                                         true};  // if true, do parabolic; if false, do linear extrapolation
        Gaudi::Property<float> m_extrapolationDistance{this, "ExtrapolationDistance", 1500.};  // default value is 1500
        Gaudi::Property<bool> m_addSectors{this, "AddSectors", false};                         // default false
        unsigned int m_ntechnologies{UINT_MAX};
        std::map<unsigned int, unsigned int> m_techToTruthNameIdx{};  // mapping the muon technology to the index of the m_truthNames vector
 
        MuonSectorMapping m_sectorMapping;
 
        mutable std::mutex m_mutex ATLAS_THREAD_SAFE;
        mutable CollectionsPerSectorVec m_collectionsPerSector ATLAS_THREAD_SAFE;
        mutable bool m_sectorSetup ATLAS_THREAD_SAFE{false};
    };
 
    struct SortHoughDataPerSector {
        bool operator()(const MuonLayerHoughTool::HoughDataPerSector* s1, const MuonLayerHoughTool::HoughDataPerSector* s2) const {
            return s2->maxEtaHits() < s1->maxEtaHits();
        }
    };
 
    inline void MuonLayerHoughTool::getSectors(const TgcClusterObj3D& tgc, std::vector<int>& sectors) const {
        getSectors(tgc.getEdge(TgcEdge::LowEtaLowPhi), sectors);
    }
 
    inline void MuonLayerHoughTool::getSectors(const Amg::Vector3D& pos, std::vector<int>& sectors) const {
        return m_sectorMapping.getSectors(pos.phi(), sectors);
    }
 
    inline double MuonLayerHoughTool::rCor(const Amg::Vector3D& pos, const Identifier& id) const {
        return m_sectorMapping.transformRToSector(pos.perp(), pos.phi(), m_idHelperSvc->sector(id));
    }
 
    inline double MuonLayerHoughTool::rCor(const MuonCluster& mm) const { return rCor(mm.globalPosition(), mm.identify()); }
 
    inline double MuonLayerHoughTool::rCor(const MdtPrepData& mm) const { return rCor(mm.globalPosition(), mm.identify()); }
 
    inline double MuonLayerHoughTool::rCor(const TgcClusterObj3D& tgc, const TgcEdge val, int sector) const {
        const Amg::Vector3D& pos{tgc.getEdge(val)};
        return m_sectorMapping.transformRToSector(pos.perp(), pos.phi(), sector);
    }
 
    inline int MuonLayerHoughTool::sublay(const Identifier& id, float /*z*/) const {
        int sublayer = 0;
        if (m_idHelperSvc->isMdt(id)) {
            sublayer = m_idHelperSvc->mdtIdHelper().tubeLayer(id) - 1;
            if (m_idHelperSvc->mdtIdHelper().multilayer(id) == 2) sublayer += 4;
        } else if (m_idHelperSvc->isMM(id)) {
            sublayer = m_idHelperSvc->mmIdHelper().gasGap(id) - 1;
            if (m_idHelperSvc->mmIdHelper().multilayer(id) == 2) sublayer += 4;
            sublayer += 600;  // type info
        } else if (m_idHelperSvc->issTgc(id)) {
            sublayer = m_idHelperSvc->stgcIdHelper().gasGap(id) - 1;
            if (m_idHelperSvc->stgcIdHelper().multilayer(id) == 2) sublayer += 4;
            sublayer += 500;  // type info
        } else if (m_idHelperSvc->isRpc(id)) {
            sublayer = m_idHelperSvc->rpcIdHelper().gasGap(id) - 1;
            if (m_idHelperSvc->rpcIdHelper().doubletR(id) == 2) sublayer += 2;
            sublayer += 100;  // type info
        } else if (m_idHelperSvc->isTgc(id)) {
            sublayer = m_idHelperSvc->tgcIdHelper().gasGap(id) - 1;
            using namespace MuonStationIndex;
            StIndex stIndex = m_idHelperSvc->stationIndex(id);
            if (stIndex == StIndex::EM) {
                // T1 gets +3; T2 gets +3+3; T3 gets +3+6; T4 gets0 (because it is also EI)
                PhiIndex phiIndex = m_idHelperSvc->phiIndex(id);
                sublayer += 3;
                if (phiIndex == PhiIndex::T2)
                    sublayer += 3;
                else if (phiIndex == PhiIndex::T3)
                    sublayer += 6;
                // float fz = fabs(z);
                // if( fz < 13350 )      sublayer += 0;
                // else if( fz < 13440 ) sublayer += 3;
                // else if( fz < 13520 ) sublayer += 6;
                // else if( fz < 13600 ) sublayer += 9;
                // else if( fz < 14650 ) sublayer += 12;
                // else if( fz < 14740 ) sublayer += 15;
                // else if( fz < 14800 ) sublayer += 18;
                // else if( fz < 14850 ) sublayer += 21;
                // else if( fz < 15070 ) sublayer += 24;
                // else if( fz < 15150 ) sublayer += 27;
                // else if( fz < 15220 ) sublayer += 30;
                // else if( fz < 15250 ) sublayer += 33;
                // else                  sublayer += 36;
            }
            sublayer += 300;  // type info
        }
        return sublayer;
    }
 
    inline double MuonLayerHoughTool::combinedPeakheight(double ph, double ph1, double ph2, double phn, double /*rot*/, int layer,
                                                         int /*region*/) const {
        if (layer == 0 && ph < 3.) return ph;
        if (layer == 1 && ph < 4.) return ph;
        if (layer == 2 && ph < 3.) return ph;
 
        if (phn > 7) ph += phn;
 
        if (layer == 0) {
            if (ph1 > 6.9) ph += 2;
            if (ph1 > 8.9) ph += 2;
            if (ph2 > 5.9) ph += 2;
            if (ph2 > 7.9) ph += 2;
        }
        if (layer == 1) {
            if (ph1 > 6.9) ph += 2;
            if (ph2 > 5.9) ph += 2;
            if (ph2 > 7.9) ph += 2;
            if (ph2 > 11.9) ph += 2;
        }
        if (layer == 2) {
            if (ph1 > 6.9) ph += 2;
            if (ph2 > 6.9) ph += 2;
            if (ph2 > 8.9) ph += 2;
        }
        return ph;
    }
 
}  // namespace Muon
#endif