diff --git a/MuonSpectrometer/MuonPhaseII/MuonDetDescr/MuonGeoModelTestR4/python/testGeoModel.py b/MuonSpectrometer/MuonPhaseII/MuonDetDescr/MuonGeoModelTestR4/python/testGeoModel.py index e9dcf9c8e9dcfffe8fb36ca5c770890ef9bb50a3..96de66cf6828af7959141bcd263a9a525283209d 100644 --- a/MuonSpectrometer/MuonPhaseII/MuonDetDescr/MuonGeoModelTestR4/python/testGeoModel.py +++ b/MuonSpectrometer/MuonPhaseII/MuonDetDescr/MuonGeoModelTestR4/python/testGeoModel.py @@ -144,6 +144,7 @@ def setupGeoR4TestCfg(args, setupSimJob = False): flags.Scheduler.ShowControlFlow = True flags.Scheduler.EnableVerboseViews = True flags.Scheduler.AutoLoadUnmetDependencies = True + flags.PerfMon.doFullMonMT = True flags.lock() diff --git a/MuonSpectrometer/MuonPhaseII/MuonPatternRecognition/MuonPatternEvent/CMakeLists.txt b/MuonSpectrometer/MuonPhaseII/MuonPatternRecognition/MuonPatternEvent/CMakeLists.txt new file mode 100644 index 0000000000000000000000000000000000000000..da4dc6630f9e07d7064fda107beb805cce00df3c --- /dev/null +++ b/MuonSpectrometer/MuonPhaseII/MuonPatternRecognition/MuonPatternEvent/CMakeLists.txt @@ -0,0 +1,10 @@ +# Copyright (C) 2002-2023 CERN for the benefit of the ATLAS collaboration + +# Declare the package name: +atlas_subdir( MuonPatternEvent ) + + +atlas_add_library( MuonPatternEvent + MuonPatternEvent/*.h Root/*.cxx + PUBLIC_HEADERS MuonPatternEvent + LINK_LIBRARIES xAODMuonPrepData Identifier MuonReadoutGeometryR4 xAODMeasurementBase MuonSpacePoint MuonStationIndexLib) diff --git a/MuonSpectrometer/MuonPhaseII/MuonPatternRecognition/MuonPatternEvent/MuonPatternEvent/HoughMaximum.h b/MuonSpectrometer/MuonPhaseII/MuonPatternRecognition/MuonPatternEvent/MuonPatternEvent/HoughMaximum.h new file mode 100644 index 0000000000000000000000000000000000000000..0d910985839f2c72e21e047b5524143c6e67b60e --- /dev/null +++ b/MuonSpectrometer/MuonPhaseII/MuonPatternRecognition/MuonPatternEvent/MuonPatternEvent/HoughMaximum.h @@ -0,0 +1,50 @@ +/* + Copyright (C) 2002-2024 CERN for the benefit of the ATLAS collaboration +*/ + +#ifndef MUONR4__HOUGHMAXIMUM__H +#define MUONR4__HOUGHMAXIMUM__H +#include <vector> +#include "xAODMeasurementBase/UncalibratedMeasurement.h" +#include "MuonSpacePoint/MuonSpacePointContainer.h" + +namespace MuonR4{ + using HoughHitType = std::shared_ptr<MuonR4::MuonSpacePoint>; + /// @brief Data class to represent one maximum in hough space. + /// For generic usage, local coordinates are referred to as "x" and "y" + /// but may correspond to any frame + class HoughMaximum{ + public: + /// @brief constructor. + /// @param x: first coordinate of maximum + /// @param y: second coordinate of maximum + /// @param counts: Size of the maximum (typically, weighted hough hit counts) - can use to sort + /// @param hits: list of hits assigned to the maximum (owned by SG). + HoughMaximum(double x, + double y, + double counts, + std::vector<HoughHitType> && hits); + /// @brief default c-tor, creates empty maximum with zero counts in the origin + HoughMaximum() = default; + /// @brief getter + /// @return the first coordinate + double getX() const; + /// @brief getter + /// @return the second coordinate + double getY() const; + /// @brief getter + /// @return the uncertainty on the first coordinate + double getCounts() const; + /// @brief getter + /// @return the hits assigned to this maximum + const std::vector<HoughHitType> & getHitsInMax() const; + + private: + double m_x{0.}; // first coordinate + double m_y{0.}; // second coordinate + double m_counts{0.}; // weighted counts + std::vector<HoughHitType> m_hitsInMax{}; // list of hits on maximum + }; +} + +#endif diff --git a/MuonSpectrometer/MuonPhaseII/MuonPatternRecognition/MuonPatternEvent/MuonPatternEvent/StationHoughMaxContainer.h b/MuonSpectrometer/MuonPhaseII/MuonPatternRecognition/MuonPatternEvent/MuonPatternEvent/StationHoughMaxContainer.h new file mode 100644 index 0000000000000000000000000000000000000000..59ea4a7cf104b3c726d6aa60c96b6e40242da0d0 --- /dev/null +++ b/MuonSpectrometer/MuonPhaseII/MuonPatternRecognition/MuonPatternEvent/MuonPatternEvent/StationHoughMaxContainer.h @@ -0,0 +1,16 @@ +/* + Copyright (C) 2002-2024 CERN for the benefit of the ATLAS collaboration +*/ + +#ifndef MUONR4_STATIONHOUGHMAXCONTAINER__H +#define MUONR4_STATIONHOUGHMAXCONTAINER__H + +#include "StationHoughMaxima.h" +#include <set> +#include "AthenaKernel/CLASS_DEF.h" +namespace MuonR4{ + using StationHoughMaxContainer = std::set<StationHoughMaxima> ; +} +CLASS_DEF( MuonR4::StationHoughMaxContainer , 1158351533 , 1 ) + +#endif diff --git a/MuonSpectrometer/MuonPhaseII/MuonPatternRecognition/MuonPatternEvent/MuonPatternEvent/StationHoughMaxima.h b/MuonSpectrometer/MuonPhaseII/MuonPatternRecognition/MuonPatternEvent/MuonPatternEvent/StationHoughMaxima.h new file mode 100644 index 0000000000000000000000000000000000000000..b8e2eb92bcf61779c4afb3fe1de55e4aa0093146 --- /dev/null +++ b/MuonSpectrometer/MuonPhaseII/MuonPatternRecognition/MuonPatternEvent/MuonPatternEvent/StationHoughMaxima.h @@ -0,0 +1,45 @@ +/* + Copyright (C) 2002-2024 CERN for the benefit of the ATLAS collaboration +*/ +#ifndef MUONR4_STATIONHOUGHMAXIMA__H +#define MUONR4_STATIONHOUGHMAXIMA__H + +#include "MuonPatternEvent/HoughMaximum.h" +#include "MuonReadoutGeometryR4/MuonChamber.h" +namespace MuonR4{ + + /// @brief Small data class to collect the hough maxima for one given station. + /// Contains a list of maxima and a station identifier. + /// Sorts by station identifier to allow set / map insertion + class StationHoughMaxima{ + public: + /// constructor + /// @param chamber: Associated chamber serving as Identifier + /// @param maxima: list of maxima (can be extended later) + StationHoughMaxima(const MuonGMR4::MuonChamber* chamber, + const std::vector<HoughMaximum> & maxima = {}); + + /// adds a maximum to the list + /// @param m: Maximum to add + void addMaximum(const HoughMaximum & m); + + /// @brief Returns the associated chamber + const MuonGMR4::MuonChamber* chamber() const; + + /// @brief getter + /// @return the maxima for this station + const std::vector<HoughMaximum> & getMaxima() const; + + /// @brief sorting operator - uses identifiers for sorting, not the maxima themselves + /// @param other: station maxima list to compare to + /// @return Comparison between the station identifiers + bool operator< (const StationHoughMaxima & other) const; + + private: + const MuonGMR4::MuonChamber* m_chamber{}; // the identifier for this station + std::vector<HoughMaximum> m_maxima{}; // the list of found maxima +}; + +} + +#endif diff --git a/MuonSpectrometer/MuonPhaseII/MuonPatternRecognition/MuonPatternEvent/Root/HoughMaximum.cxx b/MuonSpectrometer/MuonPhaseII/MuonPatternRecognition/MuonPatternEvent/Root/HoughMaximum.cxx new file mode 100644 index 0000000000000000000000000000000000000000..82686e253673561eb933966172e0e3af5218ac23 --- /dev/null +++ b/MuonSpectrometer/MuonPhaseII/MuonPatternRecognition/MuonPatternEvent/Root/HoughMaximum.cxx @@ -0,0 +1,25 @@ +/* + Copyright (C) 2002-2024 CERN for the benefit of the ATLAS collaboration +*/ +#include "MuonPatternEvent/HoughMaximum.h" + +MuonR4::HoughMaximum::HoughMaximum(double x, double y, double counts, std::vector<HoughHitType > && hits): + m_x(x), + m_y(y), + m_counts(counts), + m_hitsInMax(hits){ + +}; +double MuonR4::HoughMaximum::getX() const{ + return m_x; +} + +double MuonR4::HoughMaximum::getY() const{ + return m_y; +} +double MuonR4::HoughMaximum::getCounts() const{ + return m_counts; +} +const std::vector<MuonR4::HoughHitType> & MuonR4::HoughMaximum::getHitsInMax() const{ + return m_hitsInMax; +} diff --git a/MuonSpectrometer/MuonPhaseII/MuonPatternRecognition/MuonPatternEvent/Root/StationHoughMaxima.cxx b/MuonSpectrometer/MuonPhaseII/MuonPatternRecognition/MuonPatternEvent/Root/StationHoughMaxima.cxx new file mode 100644 index 0000000000000000000000000000000000000000..aa9d439da3ca68937160cb322d4f2ccd8409e023 --- /dev/null +++ b/MuonSpectrometer/MuonPhaseII/MuonPatternRecognition/MuonPatternEvent/Root/StationHoughMaxima.cxx @@ -0,0 +1,26 @@ +/* + Copyright (C) 2002-2024 CERN for the benefit of the ATLAS collaboration +*/ + +#include "MuonPatternEvent/StationHoughMaxima.h" +#include "MuonReadoutGeometryR4/MuonDetectorManager.h" +namespace MuonR4 { + StationHoughMaxima::StationHoughMaxima(const MuonGMR4::MuonChamber* chamber, + const std::vector<HoughMaximum> & maxima): + m_chamber{chamber}, + m_maxima(maxima){ + } + void StationHoughMaxima::addMaximum(const HoughMaximum & m){ + m_maxima.push_back(m); + } + const MuonGMR4::MuonChamber* StationHoughMaxima::chamber() const { + return m_chamber; + } + const std::vector<MuonR4::HoughMaximum> & MuonR4::StationHoughMaxima::getMaxima() const{ + return m_maxima; + } + bool MuonR4::StationHoughMaxima::operator< (const StationHoughMaxima & other) const{ + using ChamberSorter = MuonGMR4::MuonDetectorManager::ChamberSorter; + return ChamberSorter{}(m_chamber, other.m_chamber); + } +} diff --git a/MuonSpectrometer/MuonPhaseII/MuonPatternRecognition/MuonPatternRecognitionAlgs/CMakeLists.txt b/MuonSpectrometer/MuonPhaseII/MuonPatternRecognition/MuonPatternRecognitionAlgs/CMakeLists.txt new file mode 100644 index 0000000000000000000000000000000000000000..772dbe74bcd86023adf8f202b03b067e33825c0d --- /dev/null +++ b/MuonSpectrometer/MuonPhaseII/MuonPatternRecognition/MuonPatternRecognitionAlgs/CMakeLists.txt @@ -0,0 +1,17 @@ +################################################################################ +# Package: MuonPatternRecognitionAlgs +################################################################################ + +# Declare the package name: +atlas_subdir( MuonPatternRecognitionAlgs ) + +# External ACTS dependency +find_package( Acts COMPONENTS Core ) + +atlas_add_component( MuonPatternRecognitionAlgs + src/components/*.cxx src/*.cxx + LINK_LIBRARIES MuonPatternEvent ActsCore AthenaKernel StoreGateLib MuonTesterTreeLib + xAODMuonPrepData MuonSpacePoint ) + +# Install files from the package: +atlas_install_python_modules( python/*.py POST_BUILD_CMD ${ATLAS_FLAKE8} ) diff --git a/MuonSpectrometer/MuonPhaseII/MuonPatternRecognition/MuonPatternRecognitionAlgs/python/MuonHoughTransformAlgConfig.py b/MuonSpectrometer/MuonPhaseII/MuonPatternRecognition/MuonPatternRecognitionAlgs/python/MuonHoughTransformAlgConfig.py new file mode 100644 index 0000000000000000000000000000000000000000..c1054650330d2f4502f6041ce2ef42663099bf95 --- /dev/null +++ b/MuonSpectrometer/MuonPhaseII/MuonPatternRecognition/MuonPatternRecognitionAlgs/python/MuonHoughTransformAlgConfig.py @@ -0,0 +1,34 @@ +# Copyright (C) 2002-2024 CERN for the benefit of the ATLAS collaboration + +from AthenaConfiguration.ComponentAccumulator import ComponentAccumulator +from AthenaConfiguration.ComponentFactory import CompFactory + + +def MuonHoughTransformAlgCfg(flags, name = "MuonHoughTransformAlg", **kwargs): + result = ComponentAccumulator() + theAlg = CompFactory.MuonR4.MuonHoughTransformAlg(name, **kwargs) + result.addEventAlgo(theAlg, primary=True) + return result + +if __name__=="__main__": + from MuonGeoModelTestR4.testGeoModel import setupGeoR4TestCfg, SetupArgParser, executeTest,setupHistSvcCfg + parser = SetupArgParser() + parser.set_defaults(nEvents = -1) + parser.set_defaults(inputFile=["/cvmfs/atlas-nightlies.cern.ch/repo/data/data-art/MuonRecRTT/R4SimHits.pool.root"]) + + args = parser.parse_args() + flags, cfg = setupGeoR4TestCfg(args) + + from xAODMuonSimHitCnv.MuonSimHitCnvCfg import xAODSimHitToMdtMeasCnvAlgCfg + cfg.merge(setupHistSvcCfg(flags,out_file=args.outRootFile,out_stream="MuonHoughTransform")) + cfg.merge(xAODSimHitToMdtMeasCnvAlgCfg(flags)) + from MuonSpacePointFormation.SpacePointFormationConfig import MuonSpacePointMakerAlgCfg + cfg.merge(MuonSpacePointMakerAlgCfg(flags)) + cfg.merge(MuonHoughTransformAlgCfg(flags)) + + # output spam reduction + cfg.getService("AthenaHiveEventLoopMgr").EventPrintoutInterval=500 + + + executeTest(cfg, args.nEvents) + diff --git a/MuonSpectrometer/MuonPhaseII/MuonPatternRecognition/MuonPatternRecognitionAlgs/src/MuonHoughEventData.h b/MuonSpectrometer/MuonPhaseII/MuonPatternRecognition/MuonPatternRecognitionAlgs/src/MuonHoughEventData.h new file mode 100644 index 0000000000000000000000000000000000000000..bc37dab6762b71e5a653262c42575bc4a49bae80 --- /dev/null +++ b/MuonSpectrometer/MuonPhaseII/MuonPatternRecognition/MuonPatternRecognitionAlgs/src/MuonHoughEventData.h @@ -0,0 +1,64 @@ +/* + Copyright (C) 2002-2024 CERN for the benefit of the ATLAS collaboration +*/ + +#ifndef MUONR4__MUONHOUGHEVENTDATA__H +#define MUONR4__MUONHOUGHEVENTDATA__H + +#include "Acts/Seeding/HoughTransformUtils.hpp" +#include "xAODMeasurementBase/UncalibratedMeasurement.h" +#include "MuonPatternEvent/HoughMaximum.h" +#include "xAODMuonPrepData/MdtDriftCircleContainer.h" +#include "xAODMuonPrepData/RpcStripContainer.h" +#include "xAODMuonPrepData/TgcStripContainer.h" +#include "xAODMuonPrepData/MMClusterContainer.h" +#include "MuonSpacePoint/MuonSpacePointContainer.h" + + +#include <map> +#include <memory> + +// specialise our hough plane to utilise uncalibrated measurements as "identifier" type +namespace MuonR4{ + using HoughPlane = Acts::HoughTransformUtils::HoughPlane<HoughHitType> ; + using Acts::HoughTransformUtils::HoughPlaneConfig; + using ActsPeakFinderForMuon = Acts::HoughTransformUtils::PeakFinders::IslandsAroundMax<HoughHitType>; + using ActsPeakFinderForMuonCfg = Acts::HoughTransformUtils::PeakFinders::IslandsAroundMaxConfig; + class HoughMaximum; + class StationIdentifier; + template <typename peakFinder_t, typename peakFinderConfig_t> + struct MuonHoughEventData_impl{ + MuonHoughEventData_impl(const ActsGeometryContext& _gctx): + gctx{_gctx} {} + const ActsGeometryContext& gctx; + std::unique_ptr<HoughPlane> houghPlane{nullptr}; + std::unique_ptr<peakFinder_t> peakFinder{nullptr}; + + struct HoughSpaceInBucket{ + const MuonSpacePointBucket* bucket{nullptr}; + static constexpr double dummyBoundary{1.e9}; + std::pair<double, double> searchWindowZ{dummyBoundary, -dummyBoundary}; + std::pair<double, double> seachWindowTanTheta{dummyBoundary, -dummyBoundary}; + + /// Update the hough space search window + void updateHoughWindow(double zMeas, double tanThetaMeas) { + searchWindowZ.first = std::min(searchWindowZ.first, zMeas); + searchWindowZ.second = std::max(searchWindowZ.second, zMeas); + + seachWindowTanTheta.first = std::min(seachWindowTanTheta.first, tanThetaMeas); + seachWindowTanTheta.second = std::max(seachWindowTanTheta.second, tanThetaMeas); + } + }; + std::vector<HoughMaximum> maxima{}; + + std::map<const MuonGMR4::MuonChamber*, std::vector<HoughSpaceInBucket>> houghSpaces{}; + Acts::HoughTransformUtils::HoughAxisRanges currAxisRanges; + std::pair<double,double> searchSpaceZ{10000000,-100000000.}; + std::pair<double,double> searchSpaceTanTheta{100000000.,-100000.}; + }; + + // for now, we use the same (default ACTS) peak finder for both the eta- and the phi-transforms + using MuonHoughEventData = MuonHoughEventData_impl<ActsPeakFinderForMuon, ActsPeakFinderForMuonCfg>; +} + +#endif diff --git a/MuonSpectrometer/MuonPhaseII/MuonPatternRecognition/MuonPatternRecognitionAlgs/src/MuonHoughHelpers.cxx b/MuonSpectrometer/MuonPhaseII/MuonPatternRecognition/MuonPatternRecognitionAlgs/src/MuonHoughHelpers.cxx new file mode 100644 index 0000000000000000000000000000000000000000..0f0341f61ef99f2081343ac346a1f0c540ed6605 --- /dev/null +++ b/MuonSpectrometer/MuonPhaseII/MuonPatternRecognition/MuonPatternRecognitionAlgs/src/MuonHoughHelpers.cxx @@ -0,0 +1,31 @@ +/* + Copyright (C) 2002-2024 CERN for the benefit of the ATLAS collaboration +*/ + +#include "MuonHoughHelpers.h" +using namespace MuonR4; + +double HoughHelpers::Eta::houghParamMdtLeft(double tanTheta, const MuonR4::HoughHitType & DC){ + return DC->positionInChamber().y() - tanTheta * DC->positionInChamber().z() - + DC->driftRadius() * std::sqrt(1 + (tanTheta*tanTheta)); // using cos(theta) = sqrt(1/[1+tan²(theta)]) + }; + // right solution +double HoughHelpers::Eta::houghParamMdtRight(double tanTheta, const MuonR4::HoughHitType & DC){ + return DC->positionInChamber().y() - tanTheta * DC->positionInChamber().z() + + DC->driftRadius() * std::sqrt(1 + (tanTheta*tanTheta)); // using cos(theta) = sqrt(1/[1+tan²(theta)]) + }; + // solution which doesn't look at the drift circle but places a hit at the tube center +double HoughHelpers::Eta::houghParamStrip(double tanTheta, const MuonR4::HoughHitType & DC){ + return DC->positionInChamber().y() - tanTheta * DC->positionInChamber().z(); + }; + +double HoughHelpers::Eta::houghWidthMdt(double /*tanTheta*/, const MuonR4::HoughHitType & DC){ + return std::min(DC->uncertainty().x() * 3., + 1.0); // scale reported errors up to at least 1mm or 3 + // times the reported error as drift circle calib not + // fully reliable at this stage + }; +double HoughHelpers::Eta::houghWidthStrip(double /*tanTheta*/, const MuonR4::HoughHitType & DC){ + return DC->uncertainty().x(); // assign full tube radius as uncertainty + }; + diff --git a/MuonSpectrometer/MuonPhaseII/MuonPatternRecognition/MuonPatternRecognitionAlgs/src/MuonHoughHelpers.h b/MuonSpectrometer/MuonPhaseII/MuonPatternRecognition/MuonPatternRecognitionAlgs/src/MuonHoughHelpers.h new file mode 100644 index 0000000000000000000000000000000000000000..c0316a152f0eb08a073d7f61f54b2b7b8631aaeb --- /dev/null +++ b/MuonSpectrometer/MuonPhaseII/MuonPatternRecognition/MuonPatternRecognitionAlgs/src/MuonHoughHelpers.h @@ -0,0 +1,46 @@ +/* + Copyright (C) 2002-2024 CERN for the benefit of the ATLAS collaboration +*/ + +#ifndef MUONR4__MUONHOUGHHELPERS__H +#define MUONR4__MUONHOUGHHELPERS__H + +#include "MuonHoughEventData.h" +#include "Acts/Seeding/HoughTransformUtils.hpp" +#include "xAODMeasurementBase/UncalibratedMeasurement.h" +#include "xAODMuonPrepData/MdtDriftCircleContainer.h" +#include "xAODMuonPrepData/RpcStripContainer.h" + +namespace MuonR4{ + namespace HoughHelpers{ + namespace Eta{ + /// @brief left-side hough parametrisation for drift circles + /// @param tanTheta the input inclination angle + /// @param dc the drift circle + /// @return the z offset needed to touch the drift radius on the left for an inclination angle tanTheta + double houghParamMdtLeft(double tanTheta, const MuonR4::HoughHitType & dc); + /// @brief right-side hough parametrisation for drift circles + /// @param tanTheta the input inclination angle + /// @param dc the drift circle + /// @return the z offset needed to touch the drift radius on the right for an inclination angle tanTheta + double houghParamMdtRight(double tanTheta, const MuonR4::HoughHitType & dc); + /// @brief tube-level hough parametrisation for drift circles + /// @param tanTheta the input inclination angle + /// @param dc the drift circle + /// @return the z offset needed to pass through the center of the drift tube for an inclination angle tanTheta + double houghParamStrip(double tanTheta, const MuonR4::HoughHitType & dc); + /// @brief tube-level hough parametrisation for drift circles + /// @param tanTheta the input inclination angle + /// @param dc the drift circle + /// @return the z offset needed to pass through the center of the drift tube for an inclination angle tanTheta + double houghWidthMdt(double tanTheta, const MuonR4::HoughHitType & dc); + /// @brief Uncertainty parametrisation when filling the hough plane using only tube information + /// @param tanTheta: the input inclination angle (not used) + /// @param dc: the drift circle + /// @return an uncertainty corresponding to the full tube radius + double houghWidthStrip(double tanTheta, const MuonR4::HoughHitType & dc); + } + } +} + +#endif diff --git a/MuonSpectrometer/MuonPhaseII/MuonPatternRecognition/MuonPatternRecognitionAlgs/src/MuonHoughTransformAlg.cxx b/MuonSpectrometer/MuonPhaseII/MuonPatternRecognition/MuonPatternRecognitionAlgs/src/MuonHoughTransformAlg.cxx new file mode 100644 index 0000000000000000000000000000000000000000..1829f6fe806e0326a94514ba558f990761bf804b --- /dev/null +++ b/MuonSpectrometer/MuonPhaseII/MuonPatternRecognition/MuonPatternRecognitionAlgs/src/MuonHoughTransformAlg.cxx @@ -0,0 +1,154 @@ +/* + Copyright (C) 2002-2024 CERN for the benefit of the ATLAS collaboration +*/ + +#include "MuonHoughTransformAlg.h" +#include "MuonHoughHelpers.h" + +#include <StoreGate/ReadCondHandle.h> +#include <MuonReadoutGeometryR4/MuonChamber.h> + +using namespace MuonR4; + + + + +MuonHoughTransformAlg::MuonHoughTransformAlg(const std::string& name, ISvcLocator* pSvcLocator): + AthReentrantAlgorithm(name, pSvcLocator){ + } + +StatusCode MuonHoughTransformAlg::initialize(){ + ATH_CHECK(m_geoCtxKey.initialize()); + ATH_CHECK(m_spacePointKey.initialize()); + ATH_CHECK(m_maxima.initialize()); + + return StatusCode::SUCCESS; +} + +template <class ContainerType> StatusCode MuonHoughTransformAlg::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<ContainerType> readHandle{key, ctx}; + ATH_CHECK(readHandle.isPresent()); + contToPush = readHandle.cptr(); + return StatusCode::SUCCESS; +} + +StatusCode MuonHoughTransformAlg::execute(const EventContext& ctx) const { + + /// read the PRDs + const MuonSpacePointContainer* spacePoints{nullptr}; + ATH_CHECK(retrieveContainer(ctx, m_spacePointKey, spacePoints)); + + // book the output container + SG::WriteHandle<StationHoughMaxContainer> writeMaxima(m_maxima, ctx); + ATH_CHECK(writeMaxima.record(std::make_unique<StationHoughMaxContainer>())); + + SG::ReadCondHandle<ActsGeometryContext> gctxHandle{m_geoCtxKey, ctx}; + ATH_CHECK(gctxHandle.isValid()); + + MuonHoughEventData data{**gctxHandle}; + + /// pre-populate the event data - sort PRDs by station + ATH_CHECK(preProcess(data, *spacePoints)); + + /// book the hough plane + ATH_CHECK(prepareHoughPlane(data)); + /// now perform the actual HT for each station + for (auto&[station, stationHoughBuckets] : data.houghSpaces){ + ATH_CHECK(prepareStation(data,station)); + for (auto & bucket : stationHoughBuckets) { + ATH_CHECK(processBucket(data, bucket)); + } + /// write the maxima we found + writeMaxima->emplace(StationHoughMaxima(station, data.maxima)); + } + return StatusCode::SUCCESS; + +} + StatusCode MuonHoughTransformAlg::preProcess(MuonHoughEventData & data, + const MuonR4::MuonSpacePointContainer & spacePoints) const { + + ATH_MSG_DEBUG("Load "<<spacePoints.size()<<" space point buckets"); + for (const MuonR4::MuonSpacePointBucket* sp : spacePoints){ + std::vector<HoughSpaceInBucket>& buckets = data.houghSpaces[sp->front()->muonChamber()]; + buckets.push_back( + HoughSpaceInBucket{sp,} + ); + /// TODO: Update search window to fine-tune tanTheta + } + return StatusCode::SUCCESS; + } + + + StatusCode MuonHoughTransformAlg::prepareHoughPlane(MuonHoughEventData & data) const{ + HoughPlaneConfig cfg; + cfg.nBinsX = 100; + cfg.nBinsY = 100; + ActsPeakFinderForMuonCfg peakFinderCfg; + peakFinderCfg.fractionCutoff = 0.7; + peakFinderCfg.threshold = 3; + peakFinderCfg.minSpacingBetweenPeaks = {0., 30.}; + data.houghPlane = std::make_unique<HoughPlane>(cfg); + data.peakFinder = std::make_unique<ActsPeakFinderForMuon>(peakFinderCfg); + + return StatusCode::SUCCESS; + } + +StatusCode MuonHoughTransformAlg::prepareStation(MuonHoughEventData & data, const MuonGMR4::MuonChamber* ) const{ + data.maxima.clear(); + return StatusCode::SUCCESS; +} + + StatusCode MuonHoughTransformAlg::processBucket(MuonHoughEventData & data, + HoughSpaceInBucket& bucket) const{ + /// tune the search space + + double chamberCenter = 0.5 * (bucket.bucket->coveredMin() + bucket.bucket->coveredMax()); + // TODO: This will become a configurable property + double targetReso = 10.; // 1 cm target resolution in y0 + // build a symmetric window around the (geometric) chamber center so that the bin width is equivalent + // to our target resolution + // TODO: Properly pass binning instead of hardcoding nBins... + double searchStart = chamberCenter - 0.5 * 100 * targetReso; + double searchEnd = chamberCenter + 0.5 * 100 * targetReso; + // Protection for very wide buckets - if the search space does not cover all of the bucket, widen the bin size + // so that we cover everything + searchStart = std::min(searchStart, bucket.bucket->coveredMin() - 5. * targetReso); + searchEnd = std::max(searchEnd, bucket.bucket->coveredMax() + 5. * targetReso); + data.currAxisRanges = Acts::HoughTransformUtils::HoughAxisRanges{-2, 2, + searchStart, searchEnd + }; + data.houghPlane->reset(); + for (const HoughHitType & hit : *(bucket.bucket)){ + fillFromSpacePoint(data, hit); + } + auto maxima = data.peakFinder->findPeaks(*(data.houghPlane), data.currAxisRanges); + if (maxima.empty()) { + return StatusCode::SUCCESS; + } + for (const auto & max : maxima) { + /// TODO: Proper weighted hit counting... + std::vector<HoughHitType> hitList; + hitList.insert(hitList.end(), max.hitIdentifiers.begin(), max.hitIdentifiers.end()); + data.maxima.emplace_back(max.x, max.y, hitList.size(), std::move(hitList)); + } + + return StatusCode::SUCCESS; + } + void MuonHoughTransformAlg::fillFromSpacePoint(MuonHoughEventData & data, const MuonR4::HoughHitType & SP) const{ + if (SP->primaryMeasurement()->type() == xAOD::UncalibMeasType::MdtDriftCircleType){ + data.houghPlane->fill<HoughHitType>(SP, data.currAxisRanges, HoughHelpers::Eta::houghParamMdtLeft,HoughHelpers::Eta::houghWidthMdt, SP, 0, 1.0); + data.houghPlane->fill<HoughHitType>(SP, data.currAxisRanges, HoughHelpers::Eta::houghParamMdtRight,HoughHelpers::Eta::houghWidthMdt, SP, 0, 1.0); + } else{ + if (SP->measuresEta()){ + data.houghPlane->fill<HoughHitType>(SP, data.currAxisRanges, HoughHelpers::Eta::houghParamStrip,HoughHelpers::Eta::houghWidthStrip,SP, 0, 1.0); + } + } + } + diff --git a/MuonSpectrometer/MuonPhaseII/MuonPatternRecognition/MuonPatternRecognitionAlgs/src/MuonHoughTransformAlg.h b/MuonSpectrometer/MuonPhaseII/MuonPatternRecognition/MuonPatternRecognitionAlgs/src/MuonHoughTransformAlg.h new file mode 100644 index 0000000000000000000000000000000000000000..4052eb698d92c3db8dc290f40d983d82c266fac4 --- /dev/null +++ b/MuonSpectrometer/MuonPhaseII/MuonPatternRecognition/MuonPatternRecognitionAlgs/src/MuonHoughTransformAlg.h @@ -0,0 +1,69 @@ +/* + Copyright (C) 2002-2024 CERN for the benefit of the ATLAS collaboration +*/ + +#ifndef MUONR4__HOUGHTRANSFORMALG__H +#define MUONR4__HOUGHTRANSFORMALG__H + +#include "MuonPatternEvent/StationHoughMaxContainer.h" + +#include "AthenaBaseComps/AthReentrantAlgorithm.h" +#include "GaudiKernel/ToolHandle.h" +#include "StoreGate/ReadHandleKey.h" +#include "StoreGate/WriteHandleKey.h" +#include "StoreGate/ReadCondHandleKey.h" +#include "xAODMuonPrepData/MdtDriftCircleContainer.h" +#include <xAODMuonPrepData/RpcStripContainer.h> +#include <xAODMuonPrepData/TgcStripContainer.h> +#include <MuonSpacePoint/MuonSpacePointContainer.h> +#include "MuonIdHelpers/IMuonIdHelperSvc.h" +#include "MuonHoughEventData.h" + +// muon includes + + +namespace MuonR4{ + + class MuonHoughTransformAlg: public AthReentrantAlgorithm{ + public: + MuonHoughTransformAlg(const std::string& name, ISvcLocator* pSvcLocator); + virtual ~MuonHoughTransformAlg() = default; + virtual StatusCode initialize() override; + virtual StatusCode execute(const EventContext& ctx) const override; + + private: + + using HoughSpaceInBucket = MuonHoughEventData::HoughSpaceInBucket; + /// Helper method to fetch data from StoreGate. If the key is empty, a nullptr is assigned to the container ptr + /// Failure is returned in cases, of non-empty keys and failed retrieval + template <class ContainerType> StatusCode retrieveContainer(const EventContext& ctx, + const SG::ReadHandleKey<ContainerType>& key, + const ContainerType* & contToPush) const; + + + StatusCode preProcess(MuonHoughEventData & data, + const MuonSpacePointContainer & driftCircles ) const; + + StatusCode prepareHoughPlane(MuonHoughEventData & data) const; + + StatusCode prepareStation(MuonHoughEventData & data, + const MuonGMR4::MuonChamber* station) const; + StatusCode processBucket(MuonHoughEventData & data, + HoughSpaceInBucket& currentBucket) const; + + void fillFromSpacePoint(MuonHoughEventData & data, + const MuonR4::HoughHitType & SP) const; + + SG::ReadHandleKey<MuonR4::MuonSpacePointContainer> m_spacePointKey{this, "SpacePointContainer", "MuonSpacePoints"}; + + // TODO: Other technologies to join once their EDM exists + + SG::WriteHandleKey<StationHoughMaxContainer> m_maxima{this, "StationHoughMaxContainer", "MuonHoughStationMaxima"}; + + SG::ReadCondHandleKey<ActsGeometryContext> m_geoCtxKey{this, "AlignmentKey", "ActsAlignment", "cond handle key"}; + + }; +} + + +#endif diff --git a/MuonSpectrometer/MuonPhaseII/MuonPatternRecognition/MuonPatternRecognitionAlgs/src/components/MuonPatternRecognitionAlgs_entries.cxx b/MuonSpectrometer/MuonPhaseII/MuonPatternRecognition/MuonPatternRecognitionAlgs/src/components/MuonPatternRecognitionAlgs_entries.cxx new file mode 100644 index 0000000000000000000000000000000000000000..a4a980a638436419c29a974dd63b0e4f8f8df494 --- /dev/null +++ b/MuonSpectrometer/MuonPhaseII/MuonPatternRecognition/MuonPatternRecognitionAlgs/src/components/MuonPatternRecognitionAlgs_entries.cxx @@ -0,0 +1,5 @@ +/* + Copyright (C) 2002-2024 CERN for the benefit of the ATLAS collaboration +*/ +#include "../MuonHoughTransformAlg.h" +DECLARE_COMPONENT(MuonR4::MuonHoughTransformAlg) diff --git a/MuonSpectrometer/MuonPhaseII/MuonSpacePointFormation/src/MuonSpacePointMakerAlg.cxx b/MuonSpectrometer/MuonPhaseII/MuonSpacePointFormation/src/MuonSpacePointMakerAlg.cxx index 83acb8b2ab8c04d01d451a28e6891c397a0d783c..d4279010b324e5cc49f2e81d109f5b229b409555 100644 --- a/MuonSpectrometer/MuonPhaseII/MuonSpacePointFormation/src/MuonSpacePointMakerAlg.cxx +++ b/MuonSpectrometer/MuonPhaseII/MuonSpacePointFormation/src/MuonSpacePointMakerAlg.cxx @@ -213,4 +213,4 @@ void MuonSpacePointMakerAlg::distributePointsAndStore(const EventContext& ctx, } -} \ No newline at end of file +} diff --git a/MuonSpectrometer/MuonPhaseII/MuonSpacePointFormation/src/MuonSpacePointMakerAlg.h b/MuonSpectrometer/MuonPhaseII/MuonSpacePointFormation/src/MuonSpacePointMakerAlg.h index f82aa8176ad5ac46c1650f1419a7a9314380f371..b61ea62304863444939d065c0d722caaf08dd0c8 100644 --- a/MuonSpectrometer/MuonPhaseII/MuonSpacePointFormation/src/MuonSpacePointMakerAlg.h +++ b/MuonSpectrometer/MuonPhaseII/MuonSpacePointFormation/src/MuonSpacePointMakerAlg.h @@ -52,13 +52,13 @@ namespace MuonR4{ SpacePointBucketVec& splittedContainer) const; - SG::ReadHandleKey<xAOD::MdtDriftCircleContainer> m_mdtKey{this, "MdtKey", "xMdtHits", + SG::ReadHandleKey<xAOD::MdtDriftCircleContainer> m_mdtKey{this, "MdtKey", "xAODMdtCircles", "Key to the uncalibrated Drift circle measurements"}; - SG::ReadHandleKey<xAOD::RpcStripContainer> m_rpcKey{this, "RpcKey", "xRpcHits", + SG::ReadHandleKey<xAOD::RpcStripContainer> m_rpcKey{this, "RpcKey", "xRpcStrips", "Key to the uncalibrated 1D rpc hits"}; - SG::ReadHandleKey<xAOD::TgcStripContainer> m_tgcKey{this, "TgcKey", "xTgcHits", + SG::ReadHandleKey<xAOD::TgcStripContainer> m_tgcKey{this, "TgcKey", "xTgcStrips", "Key to the uncalibrated 1D tgc hits"}; SG::ReadCondHandleKey<ActsGeometryContext> m_geoCtxKey{this, "AlignmentKey", "ActsAlignment", "cond handle key"}; diff --git a/MuonSpectrometer/MuonPhaseII/MuonValidation/MuonPatternRecognitionTest/CMakeLists.txt b/MuonSpectrometer/MuonPhaseII/MuonValidation/MuonPatternRecognitionTest/CMakeLists.txt new file mode 100644 index 0000000000000000000000000000000000000000..8517c71967bc923343352a54fa78aa3b6efbc797 --- /dev/null +++ b/MuonSpectrometer/MuonPhaseII/MuonValidation/MuonPatternRecognitionTest/CMakeLists.txt @@ -0,0 +1,18 @@ +################################################################################ +# Package: MuonPatternRecognitionTest +################################################################################ + +# Declare the package name: +atlas_subdir( MuonPatternRecognitionTest ) + + +find_package( ROOT COMPONENTS Gpad Graf Core Tree MathCore Hist RIO pthread MathMore Minuit Minuit2 Matrix Physics HistPainter Rint Graf3d Html Postscript Gui GX11TTF GX11 ) + +atlas_add_component( MuonPatternRecognitionTest + src/components/*.cxx src/*.cxx + INCLUDE_DIRS ${ROOT_INCLUDE_DIRS} + LINK_LIBRARIES AthenaKernel StoreGateLib MuonTesterTreeLib + xAODMuonSimHit xAODMuonPrepData MuonPatternEvent MuonReadoutGeometryR4 ${ROOT_LIBRARIES} ) + +# Install files from the package: +atlas_install_python_modules( python/*.py POST_BUILD_CMD ${ATLAS_FLAKE8} ) diff --git a/MuonSpectrometer/MuonPhaseII/MuonValidation/MuonPatternRecognitionTest/python/MdtEtaTransformTesterConfig.py b/MuonSpectrometer/MuonPhaseII/MuonValidation/MuonPatternRecognitionTest/python/MdtEtaTransformTesterConfig.py new file mode 100644 index 0000000000000000000000000000000000000000..a21775ec8fc93e040ea4b1fa492a543136a14515 --- /dev/null +++ b/MuonSpectrometer/MuonPhaseII/MuonValidation/MuonPatternRecognitionTest/python/MdtEtaTransformTesterConfig.py @@ -0,0 +1,40 @@ +# Copyright (C) 2002-2023 CERN for the benefit of the ATLAS collaboration + +from AthenaConfiguration.ComponentAccumulator import ComponentAccumulator +from AthenaConfiguration.ComponentFactory import CompFactory + + +def MdtEtaTransformTesterCfg(flags, name = "MdtEtaTransformTester", **kwargs): + result = ComponentAccumulator() + theAlg = CompFactory.MuonValR4.MdtEtaTransformTester(name, **kwargs) + result.addEventAlgo(theAlg, primary=True) + return result + +if __name__=="__main__": + from MuonGeoModelTestR4.testGeoModel import setupGeoR4TestCfg, SetupArgParser, executeTest,setupHistSvcCfg + parser = SetupArgParser() + parser.set_defaults(nEvents = -1) + parser.set_defaults(inputFile=["/cvmfs/atlas-nightlies.cern.ch/repo/data/data-art/MuonRecRTT/R4SimHits.pool.root"]) + + args = parser.parse_args() + flags, cfg = setupGeoR4TestCfg(args) + + from xAODMuonSimHitCnv.MuonSimHitCnvCfg import MuonSimHitToMeasurementCfg + from MuonPatternRecognitionAlgs.MuonHoughTransformAlgConfig import MuonHoughTransformAlgCfg + from PerfMonComps.PerfMonCompsConfig import PerfMonMTSvcCfg + # from PerfMonVTune.PerfMonVTuneConfig import VTuneProfilerServiceCfg + cfg.merge(setupHistSvcCfg(flags,out_file=args.outRootFile,out_stream="MuonEtaHoughTransformTest")) + cfg.merge(MuonSimHitToMeasurementCfg(flags)) + from MuonSpacePointFormation.SpacePointFormationConfig import MuonSpacePointMakerAlgCfg + cfg.merge(MuonSpacePointMakerAlgCfg(flags)) + cfg.merge(MuonHoughTransformAlgCfg(flags)) + cfg.merge(MdtEtaTransformTesterCfg(flags)) + cfg.merge(PerfMonMTSvcCfg(flags)) + # cfg.merge(VTuneProfilerServiceCfg(flags, ProfiledAlgs=["MuonHoughTransformAlg"])) + + # output spam reduction + cfg.getService("AthenaHiveEventLoopMgr").EventPrintoutInterval=500 + + + executeTest(cfg, args.nEvents) + \ No newline at end of file diff --git a/MuonSpectrometer/MuonPhaseII/MuonValidation/MuonPatternRecognitionTest/src/MdtEtaTransformTester.cxx b/MuonSpectrometer/MuonPhaseII/MuonValidation/MuonPatternRecognitionTest/src/MdtEtaTransformTester.cxx new file mode 100644 index 0000000000000000000000000000000000000000..5ad96137caac4002cf199e05d9bd6ba4da5bbf9d --- /dev/null +++ b/MuonSpectrometer/MuonPhaseII/MuonValidation/MuonPatternRecognitionTest/src/MdtEtaTransformTester.cxx @@ -0,0 +1,348 @@ +/* + Copyright (C) 2002-2024 CERN for the benefit of the ATLAS collaboration +*/ + +#include "MdtEtaTransformTester.h" +#include "MuonReadoutGeometryR4/MuonChamber.h" +#include "StoreGate/ReadCondHandle.h" +#include "TCanvas.h" +#include "TLine.h" +#include "TArrow.h" +#include "TMarker.h" +#include "TEllipse.h" +#include "TLatex.h" + + + +namespace MuonValR4 { + + + + MdtEtaTransformTester::MdtEtaTransformTester(const std::string& name, + ISvcLocator* pSvcLocator) + : AthHistogramAlgorithm(name, pSvcLocator) {} + + + StatusCode MdtEtaTransformTester::initialize() { + ATH_CHECK(m_geoCtxKey.initialize()); + ATH_CHECK(m_inSimHitKeys.initialize()); + ATH_CHECK(m_inHoughMaximaKey.initialize()); + ATH_CHECK(m_spacePointKey.initialize()); + ATH_CHECK(m_tree.init(this)); + ATH_CHECK(m_idHelperSvc.retrieve()); + ATH_CHECK(detStore()->retrieve(m_r4DetMgr)); + ATH_MSG_DEBUG("Succesfully initialised"); + return StatusCode::SUCCESS; + } + + StatusCode MdtEtaTransformTester::finalize() { + ATH_CHECK(m_tree.write()); + return StatusCode::SUCCESS; + } + + Amg::Transform3D MdtEtaTransformTester::toChamberTrf(const ActsGeometryContext& gctx, + const Identifier& hitId) const { + const MuonGMR4::MuonReadoutElement* reElement = m_r4DetMgr->getReadoutElement(hitId); + //transform from local (w.r.t tube's frame) to global (ATLAS frame) and then to chamber's frame + const MuonGMR4::MuonChamber* muonChamber = reElement->getChamber(); + /// Mdt tubes have all their own transform while for the strip detectors there's one transform per layer + const IdentifierHash trfHash = reElement->detectorType() == ActsTrk::DetectorType::Mdt ? + reElement->measurementHash(hitId) : reElement->layerHash(hitId); + return muonChamber->globalToLocalTrans(gctx) * reElement->localToGlobalTrans(gctx, trfHash); + } + + StatusCode MdtEtaTransformTester::execute() { + + const EventContext & context = Gaudi::Hive::currentContext(); + SG::ReadCondHandle<ActsGeometryContext> gctxHandle{m_geoCtxKey, context}; + ATH_CHECK(gctxHandle.isValid()); + const ActsGeometryContext& gctx{**gctxHandle}; + + // retrieve the two input collections + + auto simHitCollections = m_inSimHitKeys.makeHandles(context); + + SG::ReadHandle<MuonR4::StationHoughMaxContainer> readHoughPeaks(m_inHoughMaximaKey, context); + ATH_CHECK(readHoughPeaks.isPresent()); + + ATH_MSG_DEBUG("Succesfully retrieved input collections"); + + // map the drift circles to identifiers. + // The fast digi should only generate one circle per tube. + std::map<std::pair<const MuonGMR4::MuonChamber*, HepMC::ConstGenParticlePtr>, std::vector<const xAOD::MuonSimHit*>> simHitMap{}; + for (auto & collection : simHitCollections){ + ATH_CHECK(collection.isPresent()); + for (const xAOD::MuonSimHit* simHit : *collection) { + const MuonGMR4::MuonReadoutElement* reElement = m_r4DetMgr->getReadoutElement(simHit->identify()); + const MuonGMR4::MuonChamber* id{reElement->getChamber()}; + auto genLink = simHit->genParticleLink(); + HepMC::ConstGenParticlePtr genParticle = nullptr; + if (genLink.isValid()){ + genParticle = genLink.cptr(); + } + /// skip empty truth matches for now + if (genParticle == nullptr) continue; + simHitMap[std::make_pair(id,genParticle)].push_back(simHit); + } + } + + std::map<const MuonGMR4::MuonChamber*, std::vector<MuonR4::HoughMaximum>> houghPeakMap{}; + for (const MuonR4::StationHoughMaxima & max : *readHoughPeaks){ + houghPeakMap.emplace(max.chamber(),max.getMaxima()); + } + + for (auto & [stationAndParticle, hits] : simHitMap){ + HepMC::ConstGenParticlePtr genParticlePtr = stationAndParticle.second; + const xAOD::MuonSimHit* simHit = hits.front(); + const Identifier ID = simHit->identify(); + const MuonGMR4::MuonReadoutElement* reElement = m_r4DetMgr->getReadoutElement(ID); + + const Amg::Transform3D toChamber{toChamberTrf(gctx, ID)}; + const Amg::Vector3D localPos{toChamber * xAOD::toEigen(simHit->localPosition())}; + Amg::Vector3D chamberDir = toChamber.linear() * xAOD::toEigen(simHit->localDirection()); + + /// Express the simulated hit in the center of the chamber + const std::optional<double> lambda = Amg::intersect<3>(localPos, chamberDir, Amg::Vector3D::UnitZ(), 0.); + Amg::Vector3D chamberPos = localPos + (*lambda)*chamberDir; + + m_evtNumber = context.eventID().event_number(); + m_out_stationName = reElement->stationName(); + m_out_stationEta = reElement->stationEta(); + m_out_stationPhi = reElement->stationPhi(); + /// Global coordinates + m_out_gen_Eta = genParticlePtr->momentum().eta(); + m_out_gen_Phi= genParticlePtr->momentum().phi(); + m_out_gen_Pt= genParticlePtr->momentum().pt(); + m_out_gen_nHits = hits.size(); + unsigned int nMdt{0}, nRpc{0}, nTgc{0}; + for (const xAOD::MuonSimHit* hit : hits){ + nMdt += m_idHelperSvc->isMdt(hit->identify()); + nRpc += m_idHelperSvc->isRpc(hit->identify()); + nTgc += m_idHelperSvc->isTgc(hit->identify()); + } + m_out_gen_nRPCHits = nRpc; + m_out_gen_nMDTHits = nMdt; + m_out_gen_nTGCHits = nTgc; + + m_out_gen_tantheta = (std::abs(chamberDir.z()) > 1.e-8 ? chamberDir.y()/chamberDir.z() : 1.e10); + m_out_gen_z0 = chamberPos.y(); + + const std::vector<MuonR4::HoughMaximum>& houghMaxima = houghPeakMap[reElement->getChamber()]; + if (houghMaxima.empty()){ + if (!m_tree.fill(context)) return StatusCode::FAILURE; + continue; + } + const MuonR4::HoughMaximum* foundMax = nullptr; + // find the best hough maximum + size_t max_hits{0}; + for (const MuonR4::HoughMaximum & max : houghMaxima){ + size_t nFound{0}; + for (const xAOD::MuonSimHit* simHit : hits) { + + for (const MuonR4::HoughHitType & hitOnMax : max.getHitsInMax()) { + if(m_idHelperSvc->isMdt(simHit->identify()) && + hitOnMax->identify() == simHit->identify()){ + ++nFound; + break; + } + //// Sim hits are expressed w.r.t to the gas gap Id. Check whether + /// the hit is in the same gas gap + else if (m_idHelperSvc->gasGapId(hitOnMax->identify()) == simHit->identify()) { + ++nFound; + } + } + } + if (nFound > max_hits){ + max_hits = nFound; + foundMax = &max; + } + } + /// Maximum could be associated to the hit + if (foundMax != nullptr){ + m_out_hasMax = true; + m_out_max_tantheta = foundMax->getX(); + m_out_max_z0 = foundMax->getY(); + m_out_max_nHits = max_hits; + unsigned int nMdt{0}, nRpc{0}, nTgc{0}; + for (const MuonR4::HoughHitType & houghSP: foundMax->getHitsInMax()){ + /// Skip all space points that don' contain any phi measurement + if (!houghSP->measuresEta()) continue; + + const xAOD::UncalibratedMeasurement* meas = houghSP->primaryMeasurement(); + switch (meas->type()) { + case xAOD::UncalibMeasType::MdtDriftCircleType: + m_max_driftCircleId.push_back(houghSP->identify()); + m_max_driftCircleTubePos.push_back(houghSP->positionInChamber()); + m_max_driftCirclRadius.push_back(houghSP->driftRadius()); + m_max_driftCircleDriftUncert.push_back(houghSP->uncertainty()[0]); + m_max_driftCircleTubeLength.push_back(houghSP->uncertainty()[1]); + ++nMdt; + break; + case xAOD::UncalibMeasType::RpcStripType: + m_max_rpcHitId.push_back(houghSP->identify()); + m_max_rpcHitPos.push_back(houghSP->positionInChamber()); + m_max_rpcHitHasPhiMeas.push_back(houghSP->measuresPhi()); + m_max_rpcHitErrorX.push_back(houghSP->uncertainty()[0]); + m_max_rpcHitErrorY.push_back(houghSP->uncertainty()[1]); + ++nRpc; + break; + case xAOD::UncalibMeasType::TgcStripType: + m_max_tgcHitId.push_back(houghSP->identify()); + m_max_tgcHitPos.push_back(houghSP->positionInChamber()); + m_max_tgcHitHasPhiMeas.push_back(houghSP->measuresPhi()); + m_max_tgcHitErrorX.push_back(houghSP->uncertainty()[0]); + m_max_tgcHitErrorY.push_back(houghSP->uncertainty()[1]); + ++nTgc; + break; + default: + ATH_MSG_WARNING("Technology "<<m_idHelperSvc->toString(houghSP->identify()) + <<" not yet implemented"); + } + } + m_out_max_nMdt = nMdt; + m_out_max_nRpc = nRpc; + m_out_max_nTgc = nTgc; + } + + if (!m_tree.fill(context)) return StatusCode::FAILURE; + } + return StatusCode::SUCCESS; + } + StatusCode MdtEtaTransformTester::drawEventDisplay(const EventContext& ctx, + const std::vector<const xAOD::MuonSimHit*>& simHits, + const MuonR4::HoughMaximum* foundMax) const { + + if (simHits.size() < 4) return StatusCode::SUCCESS; + + SG::ReadCondHandle<ActsGeometryContext> gctxHandle{m_geoCtxKey, ctx}; + ATH_CHECK(gctxHandle.isValid()); + const ActsGeometryContext& gctx{**gctxHandle}; + + auto readSpacePoints = SG::makeHandle(m_spacePointKey, ctx); + ATH_CHECK(readSpacePoints.isPresent()); + + double zmin{1.e9}, zmax{-1.e9}, ymin{1.e9}, ymax{-1.e9}; + std::vector<std::pair<double,double>> shPos{}, shDir{}; + + const MuonGMR4::MuonChamber* refChamber = m_r4DetMgr->getChamber(simHits[0]->identify()); + + for (const xAOD::MuonSimHit* thisHit: simHits){ + const Identifier thisID = thisHit->identify(); + const Amg::Transform3D toChamber = toChamberTrf(gctx, thisID); + + const Amg::Vector3D localPos{toChamber * xAOD::toEigen(thisHit->localPosition())}; + const Amg::Vector3D chamberDir = toChamber.linear() * xAOD::toEigen(thisHit->localDirection()); + + shPos.push_back(std::make_pair(localPos.y(), localPos.z())); + shDir.push_back(std::make_pair(chamberDir.y(), chamberDir.z())); + zmin = std::min(localPos.z(), zmin); + zmax = std::max(localPos.z(), zmax); + ymin = std::min(localPos.y(), ymin); + ymax = std::max(localPos.y(), ymax); + } + + TCanvas myCanvas("can","can",800,600); + myCanvas.cd(); + + double width = (ymax - ymin)*1.1; + double height = (zmax - zmin)*1.1; + if (height > width) width = height; + else height = width; + + double y0 = 0.5 * (ymax + ymin) - 0.5 * width; + double z0 = 0.5 * (zmax + zmin) - 0.5 * height; + double y1 = 0.5 * (ymax + ymin) + 0.5 * width; + double z1 = 0.5 * (zmax + zmin) + 0.5 * height; + auto frame = myCanvas.DrawFrame(y0,z0,y1,z1); + double frameWidth = frame->GetXaxis()->GetXmax() - frame->GetXaxis()->GetXmin(); + + std::vector<std::unique_ptr<TObject>> primitives; + for (size_t h = 0; h < shPos.size(); ++h){ + auto l = std::make_unique<TArrow>(shPos.at(h).first, shPos.at(h).second,shPos.at(h).first + 0.3 * frameWidth * shDir.at(h).first, shPos.at(h).second + 0.3 * frameWidth * shDir.at(h).second); + l->SetLineStyle(kDotted); + l->Draw(); + primitives.emplace_back(std::move(l)); + auto m = std::make_unique<TMarker>(shPos.at(h).first, shPos.at(h).second,kFullDotLarge); + m->Draw(); + primitives.emplace_back(std::move(m)); + } + + for (auto spbucket : *readSpacePoints) { + if (spbucket->muonChamber() != refChamber) continue; + std::set<MuonR4::HoughHitType> hitsOnMax{}; + if (foundMax){ + hitsOnMax.insert(foundMax->getHitsInMax().begin(), foundMax->getHitsInMax().end()); + } + for (auto & hit : *spbucket){ + ATH_MSG_DEBUG( " HIT @ "<<Amg::toString(hit->positionInChamber())<<" "<<m_idHelperSvc->toString(hit->identify())<<" with r = "<<hit->driftRadius()); + /// Space point is a mdt space point + if (hit->driftRadius() > 1e-6) { + auto ell = std::make_unique<TEllipse>(hit->positionInChamber().y(), hit->positionInChamber().z(),hit->driftRadius()); + ell->SetLineColor(kRed); + ell->SetFillColor(kRed); + if (hitsOnMax.count(hit)) { + ell->SetFillStyle(1001); + ell->SetFillColorAlpha(ell->GetFillColor(), 0.8); + }else{ + ell->SetFillStyle(0); + } + ell->Draw(); + primitives.emplace_back(std::move(ell)); + } else { + auto m = std::make_unique<TEllipse>(hit->positionInChamber().y(), hit->positionInChamber().z(), hit->uncertainty().x(), hit->uncertainty().x()); + m->SetLineColor(kBlue); + m->SetFillColor(kBlue); + m->SetFillStyle(0); + + if (hit->measuresPhi() && hit->measuresEta()) { + m->SetLineColor(kViolet); + m->SetFillColor(kViolet); + } else if (hit->measuresPhi()) { + m->SetLineColor(kGreen); + m->SetFillColor(kGreen); + } + /// Fill the ellipse if the hit is on the hough maximum + if (hitsOnMax.count(hit)) { + m->SetFillStyle(1001); + m->SetFillColorAlpha(m->GetFillColor(), 0.8); + } + m->Draw(); + primitives.emplace_back(std::move(m)); + } + } + } + + std::stringstream legendLabel{}; + legendLabel<<"Evt "<<ctx.evt()<<" station: "<<m_idHelperSvc->mdtIdHelper().stationNameString(refChamber->stationName()); + legendLabel<<"eta: "<<refChamber->stationEta()<<", phi: "<<refChamber->stationPhi(); + legendLabel<<", found maximum: "<<( foundMax ? "si" : "no"); + + + TLatex tl( 0.15,0.8,legendLabel.str().c_str()); + tl.SetNDC(); + tl.SetTextFont(53); + tl.SetTextSize(18); + tl.Draw(); + if (foundMax) { + auto mrk = std::make_unique<TMarker>( foundMax->getY(), 0., kFullTriangleUp); + mrk->SetMarkerSize(1); + mrk->SetMarkerColor(kOrange-3); + mrk->Draw(); + primitives.emplace_back(std::move(mrk)); + auto trajectory = std::make_unique<TArrow>( foundMax->getY(), 0., foundMax->getY() + 0.3 * frameWidth * m_out_max_tantheta.getVariable(), 0.3 * frameWidth); + trajectory->SetLineColor(kOrange-3); + trajectory->Draw(); + primitives.push_back(std::move(trajectory)); + } + static std::atomic<unsigned int> pdfCounter{0}; + std::stringstream pdfName{}; + pdfName<<"HoughEvt_"<<ctx.evt()<<"_"<<(++pdfCounter) + <<m_idHelperSvc->mdtIdHelper().stationNameString(refChamber->stationName())<<"_" + <<refChamber->stationEta()<<"_"<<refChamber->stationPhi()<<".pdf"; + + myCanvas.SaveAs(pdfName.str().c_str()); + primitives.clear(); + return StatusCode::SUCCESS; + } + +} // namespace MuonValR4 diff --git a/MuonSpectrometer/MuonPhaseII/MuonValidation/MuonPatternRecognitionTest/src/MdtEtaTransformTester.h b/MuonSpectrometer/MuonPhaseII/MuonValidation/MuonPatternRecognitionTest/src/MdtEtaTransformTester.h new file mode 100644 index 0000000000000000000000000000000000000000..8bc1c104d5d1b52b07c3b02d0323d180449cb82f --- /dev/null +++ b/MuonSpectrometer/MuonPhaseII/MuonValidation/MuonPatternRecognitionTest/src/MdtEtaTransformTester.h @@ -0,0 +1,110 @@ +/* + Copyright (C) 2002-2024 CERN for the benefit of the ATLAS collaboration +*/ + +#ifndef MUONFASTDIGITEST_MUONVALR4_MdtEtaTransformTester_H +#define MUONFASTDIGITEST_MUONVALR4_MdtEtaTransformTester_H + +// Framework includes +#include "AthenaBaseComps/AthHistogramAlgorithm.h" +#include "StoreGate/ReadHandleKey.h" +#include "StoreGate/ReadHandleKeyArray.h" +#include "StoreGate/ReadCondHandleKey.h" + +// EDM includes +#include "xAODMuonSimHit/MuonSimHitContainer.h" +#include "xAODMuonPrepData/MdtDriftCircleContainer.h" +#include "xAODMuonPrepData/RpcStripContainer.h" +#include "xAODMuonPrepData/TgcStripContainer.h" +#include "MuonPatternEvent/StationHoughMaxContainer.h" +#include <MuonReadoutGeometryR4/MuonDetectorManager.h> + +// muon includes +#include "MuonIdHelpers/IMuonIdHelperSvc.h" +#include "MuonTesterTree/MuonTesterTree.h" +#include "MuonTesterTree/ThreeVectorBranch.h" +#include "MuonTesterTree/IdentifierBranch.h" + + +/// @brief Lightweight algorithm to read xAOD MDT sim hits and +/// (fast-digitised) drift circles from SG and fill a +/// validation NTuple with identifier and drift circle info. + +namespace MuonValR4{ + + class MdtEtaTransformTester : public AthHistogramAlgorithm { + public: + MdtEtaTransformTester(const std::string& name, ISvcLocator* pSvcLocator); + virtual ~MdtEtaTransformTester() = default; + + virtual StatusCode initialize() override; + virtual StatusCode execute() override; + virtual StatusCode finalize() override; + + private: + Amg::Transform3D toChamberTrf(const ActsGeometryContext& gctx, + const Identifier& hitId) const; + StatusCode drawEventDisplay(const EventContext& ctx, + const std::vector<const xAOD::MuonSimHit*>& simHits, + const MuonR4::HoughMaximum* foundMax) const; + + // MDT sim hits in xAOD format + SG::ReadHandleKeyArray<xAOD::MuonSimHitContainer> m_inSimHitKeys {this, "SimHitKeys",{ "xMdtSimHits","xRpcSimHits","xTgcSimHits"}, "xAOD SimHit collections"}; + + SG::ReadHandleKey<MuonR4::StationHoughMaxContainer> m_inHoughMaximaKey{this, "StationHoughMaxContainer", "MuonHoughStationMaxima"}; + SG::ReadHandleKey<MuonR4::MuonSpacePointContainer> m_spacePointKey{this, "SpacePointContainer", "MuonSpacePoints"}; + + SG::ReadCondHandleKey<ActsGeometryContext> m_geoCtxKey{this, "AlignmentKey", "ActsAlignment", "cond handle key"}; + + ServiceHandle<Muon::IMuonIdHelperSvc> m_idHelperSvc{this, "MuonIdHelperSvc", "Muon::MuonIdHelperSvc/MuonIdHelperSvc"}; + const MuonGMR4::MuonDetectorManager* m_r4DetMgr{nullptr}; + + // // output tree - allows to compare the sim and fast-digitised hits + MuonVal::MuonTesterTree m_tree{"MuonEtaHoughTest","MuonEtaHoughTransformTest"}; + MuonVal::ScalarBranch<Long64_t>& m_evtNumber{m_tree.newScalar<Long64_t>("eventNumber")}; + MuonVal::ScalarBranch<int>& m_out_stationName{m_tree.newScalar<int>("stationName")}; + MuonVal::ScalarBranch<int>& m_out_stationEta{m_tree.newScalar<int>("stationEta")}; + MuonVal::ScalarBranch<int>& m_out_stationPhi{m_tree.newScalar<int>("stationPhi")}; + MuonVal::ScalarBranch<float>& m_out_gen_Eta{m_tree.newScalar<float>("genEta")}; + MuonVal::ScalarBranch<float>& m_out_gen_Phi{m_tree.newScalar<float>("genPhi")}; + MuonVal::ScalarBranch<float>& m_out_gen_Pt{m_tree.newScalar<float>("genPt")}; + MuonVal::ScalarBranch<unsigned int>& m_out_gen_nHits{m_tree.newScalar<unsigned int>("genNHits")}; + MuonVal::ScalarBranch<unsigned int>& m_out_gen_nRPCHits{m_tree.newScalar<unsigned int>("genNRpcHits",0)}; + MuonVal::ScalarBranch<unsigned int>& m_out_gen_nMDTHits{m_tree.newScalar<unsigned int>("genNMdtHits",0)}; + MuonVal::ScalarBranch<unsigned int>& m_out_gen_nTGCHits{m_tree.newScalar<unsigned int>("genNTgcHits",0)}; + + MuonVal::ScalarBranch<float>& m_out_gen_tantheta{m_tree.newScalar<float>("genTanTheta", 0.0)}; + MuonVal::ScalarBranch<float>& m_out_gen_z0{m_tree.newScalar<float>("genZ0", 0.0)}; + MuonVal::ScalarBranch<bool>& m_out_hasMax {m_tree.newScalar<bool>("hasMax", false)}; + MuonVal::ScalarBranch<float>& m_out_max_tantheta{m_tree.newScalar<float>("maxTanTheta", 0.0)}; + MuonVal::ScalarBranch<float>& m_out_max_z0{m_tree.newScalar<float>("maxZ0", 0.0)}; + + MuonVal::ScalarBranch<unsigned int>& m_out_max_nHits{m_tree.newScalar<unsigned int>("maxNHits", 0)}; + MuonVal::ScalarBranch<unsigned int>& m_out_max_nMdt{m_tree.newScalar<unsigned int>("maxNMdtHits", 0)}; + MuonVal::ScalarBranch<unsigned int>& m_out_max_nRpc{m_tree.newScalar<unsigned int>("maxNRpcHits", 0)}; + MuonVal::ScalarBranch<unsigned int>& m_out_max_nTgc{m_tree.newScalar<unsigned int>("maxNTgcHits", 0)}; + /// Dump of the Mdt hits on maximum + MuonVal::MdtIdentifierBranch m_max_driftCircleId{m_tree, "maxMdtId"}; + MuonVal::VectorBranch<float>& m_max_driftCirclRadius{m_tree.newVector<float>("maxMdtDriftR")}; + MuonVal::ThreeVectorBranch m_max_driftCircleTubePos{m_tree,"maxMdtTubePos"}; + MuonVal::VectorBranch<float>& m_max_driftCircleDriftUncert{m_tree.newVector<float>("maxMdtUncertDriftR")}; + MuonVal::VectorBranch<float>& m_max_driftCircleTubeLength{m_tree.newVector<float>("maxMdtUncertWire")}; + + + /// @brief Branches to access the space points in the maximum + MuonVal::RpcIdentifierBranch m_max_rpcHitId{m_tree, "maxRpcId"}; + MuonVal::ThreeVectorBranch m_max_rpcHitPos{m_tree,"maxRpcHitPos"}; + MuonVal::VectorBranch<bool> & m_max_rpcHitHasPhiMeas{m_tree.newVector<bool>("maxRpcHasPhiMeas")}; + MuonVal::VectorBranch<float>& m_max_rpcHitErrorX{m_tree.newVector<float>("maxRpcEtaMeasError")}; + MuonVal::VectorBranch<float>& m_max_rpcHitErrorY{m_tree.newVector<float>("maxRpcPhiMeasError")}; + + MuonVal::TgcIdentifierBranch m_max_tgcHitId{m_tree, "maxTgcId"}; + MuonVal::ThreeVectorBranch m_max_tgcHitPos{m_tree,"maxTgcHitPos"}; + MuonVal::VectorBranch<bool> & m_max_tgcHitHasPhiMeas{m_tree.newVector<bool>("maxTgcHasPhiMeas")}; + MuonVal::VectorBranch<float>& m_max_tgcHitErrorX{m_tree.newVector<float>("maxTgcEtaMeasError")}; + MuonVal::VectorBranch<float>& m_max_tgcHitErrorY{m_tree.newVector<float>("maxTgcPhiMeasError")}; + + }; +} + +#endif // MUONFASTDIGITEST_MUONVALR4_MdtEtaTransformTester_H diff --git a/MuonSpectrometer/MuonPhaseII/MuonValidation/MuonPatternRecognitionTest/src/components/MuonPatternRecognitionTest_entries.cxx b/MuonSpectrometer/MuonPhaseII/MuonValidation/MuonPatternRecognitionTest/src/components/MuonPatternRecognitionTest_entries.cxx new file mode 100644 index 0000000000000000000000000000000000000000..601458c6c0752e40098c88d9e424f496b9551a39 --- /dev/null +++ b/MuonSpectrometer/MuonPhaseII/MuonValidation/MuonPatternRecognitionTest/src/components/MuonPatternRecognitionTest_entries.cxx @@ -0,0 +1,7 @@ +/* + Copyright (C) 2002-2024 CERN for the benefit of the ATLAS collaboration +*/ + +#include "../MdtEtaTransformTester.h" + +DECLARE_COMPONENT(MuonValR4::MdtEtaTransformTester)