TrackFinding/src/MuonHoughSeeder.cpp

0001 // This file is part of the ACTS project.
0002 //
0003 // Copyright (C) 2016 CERN for the benefit of the ACTS project
0004 //
0005 // This Source Code Form is subject to the terms of the Mozilla Public
0006 // License, v. 2.0. If a copy of the MPL was not distributed with this
0007 // file, You can obtain one at https://mozilla.org/MPL/2.0/.
0008
0009 #include "ActsExamples/TrackFinding/MuonHoughSeeder.hpp"
0010
0011 #include "Acts/Definitions/Units.hpp"
0012 #include "ActsExamples/EventData/MuonHoughMaximum.hpp"
0013 #include "ActsExamples/EventData/MuonSegment.hpp"
0014 #include "ActsExamples/EventData/MuonSpacePoint.hpp"
0015
0016 #include <algorithm>
0017 #include <cmath>
0018 #include <format>
0019 #include <stdexcept>
0020
0021 namespace ActsExamples {
0022
0023 // left solution
0024 auto etaHoughParamDC_left = [](double tanAlpha, const MuonSpacePoint& DC) {
0025   return DC.localPosition().y() - tanAlpha * DC.localPosition().z() -
0026          DC.driftRadius() * std::sqrt(1. + tanAlpha * tanAlpha);
0027 };
0028 // right solution
0029 auto etaHoughParamDC_right = [](double tanAlpha, const MuonSpacePoint& DC) {
0030   return DC.localPosition().y() - tanAlpha * DC.localPosition().z() +
0031          DC.driftRadius() * std::sqrt(1. + tanAlpha * tanAlpha);
0032 };
0033
0034 // create the function parametrising the drift radius uncertainty
0035 auto houghWidth_fromDC = [](double, const MuonSpacePoint& DC) {
0036   // scale reported errors up to at least 1mm or 3 times the reported error as
0037   // drift circle calib not fully reliable at this stage
0038   return std::min(std::sqrt(DC.covariance()[Acts::eY]) * 3., 1.0);
0039 };
0040
0041 /// strip solution
0042 auto etaHoughParam_strip = [](double tanAlpha, const MuonSpacePoint& strip) {
0043   return strip.localPosition().y() - tanAlpha * strip.localPosition().z();
0044 };
0045 auto phiHoughParam_strip = [](double tanBeta, const MuonSpacePoint& strip) {
0046   return strip.localPosition().x() - tanBeta * strip.localPosition().z();
0047 };
0048
0049 /// @brief Strip uncertainty
0050 auto etaHoughWidth_strip = [](double, const MuonSpacePoint& strip) {
0051   return std::sqrt(strip.covariance()[Acts::eY]) * 3.;
0052 };
0053 auto phiHoughWidth_strip = [](double, const MuonSpacePoint& strip) {
0054   return std::sqrt(strip.covariance()[Acts::eX]) * 3.;
0055 };
0056
0057 MuonHoughSeeder::MuonHoughSeeder(const Config& cfg,
0058                                  std::unique_ptr<const Acts::Logger> logger)
0059     : IAlgorithm("MuonHoughSeeder", std::move(logger)), m_cfg(cfg) {
0060   if (m_cfg.inSpacePoints.empty()) {
0061     throw std::invalid_argument(
0062         "MuonHoughSeeder: Missing drift circle collection");
0063   }
0064   if (m_cfg.inTruthSegments.empty()) {
0065     throw std::invalid_argument(
0066         "MuonHoughSeeder: Missing truth segment collection");
0067   }
0068
0069   m_inputSpacePoints.initialize(m_cfg.inSpacePoints);
0070   m_inputTruthSegs.initialize(m_cfg.inTruthSegments);
0071   m_outputMaxima.initialize(m_cfg.outHoughMax);
0072 }
0073
0074 MuonHoughSeeder::~MuonHoughSeeder() = default;
0075
0076 ProcessCode MuonHoughSeeder::execute(const AlgorithmContext& ctx) const {
0077   // read the hits and circles
0078   const MuonSpacePointContainer& gotSpacePoints = m_inputSpacePoints(ctx);
0079
0080   // configure the binning of the hough plane
0081   Acts::HoughTransformUtils::HoughPlaneConfig etaPlaneCfg;
0082   etaPlaneCfg.nBinsX = m_cfg.nBinsTanTheta;
0083   etaPlaneCfg.nBinsY = m_cfg.nBinsY0;
0084
0085   Acts::HoughTransformUtils::HoughPlaneConfig phiPlaneCfg;
0086   phiPlaneCfg.nBinsX = m_cfg.nBinsTanPhi;
0087   phiPlaneCfg.nBinsY = m_cfg.nBinsX0;
0088
0089   // instantiate the hough plane
0090   HoughPlane_t etaPlane{etaPlaneCfg};
0091   HoughPlane_t phiPlane{phiPlaneCfg};
0092   MuonHoughMaxContainer outMaxima{};
0093
0094   for (const MuonSpacePointContainer::value_type& bucket : gotSpacePoints) {
0095     MuonHoughMaxContainer etaMaxima = constructEtaMaxima(ctx, bucket, etaPlane);
0096     ACTS_VERBOSE(__func__ << "() - " << __LINE__ << " Found "
0097                           << etaMaxima.size() << " eta maxima");
0098     MuonHoughMaxContainer stationMaxima =
0099         extendMaximaWithPhi(ctx, std::move(etaMaxima), phiPlane);
0100     ACTS_VERBOSE(__func__ << "() - " << __LINE__ << " Extended the maxima to "
0101                           << stationMaxima.size() << " station maxima");
0102     outMaxima.insert(outMaxima.end(),
0103                      std::make_move_iterator(stationMaxima.begin()),
0104                      std::make_move_iterator(stationMaxima.end()));
0105   }
0106
0107   m_outputMaxima(ctx, std::move(outMaxima));
0108
0109   return ProcessCode::SUCCESS;
0110 }
0111
0112 MuonHoughMaxContainer MuonHoughSeeder::constructEtaMaxima(
0113     const AlgorithmContext& ctx, const MuonSpacePointBucket& bucket,
0114     HoughPlane_t& plane) const {
0115   MuonHoughMaxContainer etaMaxima{};
0116   AxisRange_t axisRanges{-3., 3., 100. * Acts::UnitConstants::m,
0117                          -100. * Acts::UnitConstants::m};
0118   /** brief Adapt the intercept axis ranges */
0119   for (const MuonSpacePoint& sp : bucket) {
0120     /** Require that the space point measures eta */
0121     if (!sp.id().measuresEta()) {
0122       continue;
0123     }
0124     const double y = sp.localPosition().y();
0125
0126     axisRanges.yMin = std::min(axisRanges.yMin, y - m_cfg.etaPlaneMarginIcept);
0127     axisRanges.yMax = std::min(axisRanges.yMax, y + m_cfg.etaPlaneMarginIcept);
0128   }
0129   /** Ranges are adapted. Now fill the hough plane */
0130   plane.reset();
0131   for (const MuonSpacePoint& sp : bucket) {
0132     /** Require that the space point measures eta */
0133     if (!sp.id().measuresEta()) {
0134       continue;
0135     }
0136     if (sp.id().technology() == MuonSpacePoint::MuonId::TechField::Mdt) {
0137       plane.fill<MuonSpacePoint>(sp, axisRanges, etaHoughParamDC_left,
0138                                  houghWidth_fromDC, &sp, sp.id().detLayer());
0139       plane.fill<MuonSpacePoint>(sp, axisRanges, etaHoughParamDC_right,
0140                                  houghWidth_fromDC, &sp, sp.id().detLayer());
0141     } else {
0142       plane.fill<MuonSpacePoint>(sp, axisRanges, etaHoughParam_strip,
0143                                  etaHoughWidth_strip, &sp, sp.id().detLayer());
0144     }
0145   }
0146   /** Extract the maxima from the peak */
0147   PeakFinderCfg_t peakFinderCfg{};
0148   peakFinderCfg.fractionCutoff = 0.7f;
0149   peakFinderCfg.threshold = 3.f;
0150   peakFinderCfg.minSpacingBetweenPeaks = {0.f, 30.f};
0151
0152   PeakFinder_t peakFinder{peakFinderCfg};
0153
0154   /// Fetch the first set of peaks
0155   const MaximumVec_t maxima = peakFinder.findPeaks(plane, axisRanges);
0156   if (maxima.empty()) {
0157     ACTS_VERBOSE(__func__ << "() - " << __LINE__
0158                           << " No maxima found for bucket");
0159     return etaMaxima;
0160   }
0161   /** Create a measurement vector of the pure phi hits to be copied on all eta
0162    * maxima */
0163   MuonHoughMaximum::HitVec phiHits{};
0164   for (const MuonSpacePoint& sp : bucket) {
0165     if (!sp.id().measuresEta()) {
0166       phiHits.emplace_back(&sp);
0167     }
0168   }
0169   for (const Maximum_t& peakMax : maxima) {
0170     MuonHoughMaximum::HitVec hits{peakMax.hitIdentifiers.begin(),
0171                                   peakMax.hitIdentifiers.end()};
0172     hits.insert(hits.end(), phiHits.begin(), phiHits.end());
0173     etaMaxima.emplace_back(peakMax.x, peakMax.y, std::move(hits));
0174   }
0175   MuonId bucketId{bucket.front().id()};
0176   bucketId.setCoordFlags(true, false);
0177   displayMaxima(ctx, bucketId, maxima, plane, axisRanges);
0178   return etaMaxima;
0179 }
0180 MuonHoughMaxContainer MuonHoughSeeder::extendMaximaWithPhi(
0181     const AlgorithmContext& ctx, MuonHoughMaxContainer&& etaMaxima,
0182     HoughPlane_t& plane) const {
0183   MuonHoughMaxContainer outMaxima{};
0184
0185   PeakFinderCfg_t peakFinderCfg{};
0186   peakFinderCfg.fractionCutoff = 0.7;
0187   peakFinderCfg.threshold = 2.;
0188   peakFinderCfg.minSpacingBetweenPeaks = {0., 30.};
0189
0190   PeakFinder_t peakFinder{peakFinderCfg};
0191
0192   for (MuonHoughMaximum& etaMax : etaMaxima) {
0193     plane.reset();
0194     AxisRange_t axisRanges{-3., 3., 100. * Acts::UnitConstants::m,
0195                            -100. * Acts::UnitConstants::m};
0196
0197     unsigned nPhiHits{0};
0198     for (const MuonSpacePoint* sp : etaMax.hits()) {
0199       if (!sp->id().measuresPhi()) {
0200         continue;
0201       }
0202       ++nPhiHits;
0203       const double x = sp->localPosition().x();
0204       axisRanges.yMax =
0205           std::min(axisRanges.yMax, x + m_cfg.phiPlaneMarginIcept);
0206       axisRanges.yMin =
0207           std::min(axisRanges.yMin, x - m_cfg.phiPlaneMarginIcept);
0208     }
0209     if (nPhiHits < 2) {
0210       outMaxima.emplace_back(std::move(etaMax));
0211       continue;
0212     }
0213     MuonHoughMaximum::HitVec etaHits{};
0214     etaHits.reserve(etaMax.hits().size());
0215     for (const MuonSpacePoint* sp : etaMax.hits()) {
0216       if (!sp->id().measuresPhi()) {
0217         etaHits.push_back(sp);
0218         continue;
0219       }
0220       plane.fill<MuonSpacePoint>(*sp, axisRanges, phiHoughParam_strip,
0221                                  phiHoughWidth_strip, sp, sp->id().detLayer());
0222     }
0223     const MaximumVec_t phiMaxima = peakFinder.findPeaks(plane, axisRanges);
0224     for (const Maximum_t& max : phiMaxima) {
0225       MuonHoughMaximum::HitVec hits{max.hitIdentifiers.begin(),
0226                                     max.hitIdentifiers.end()};
0227       hits.insert(hits.end(), etaHits.begin(), etaHits.end());
0228       outMaxima.emplace_back(max.x, max.y, etaMax.tanBeta(),
0229                              etaMax.interceptY(), std::move(hits));
0230     }
0231     MuonId bucketId{etaMax.hits().front()->id()};
0232     bucketId.setCoordFlags(false, true);
0233     displayMaxima(ctx, bucketId, phiMaxima, plane, axisRanges);
0234   }
0235   return outMaxima;
0236 }
0237
0238 void MuonHoughSeeder::displayMaxima(const AlgorithmContext& ctx,
0239                                     const MuonId& bucketId,
0240                                     const MaximumVec_t& maxima,
0241                                     const HoughPlane_t& plane,
0242                                     const AxisRange_t& axis) const {
0243   if (!m_cfg.dumpVisualization || !m_cfg.visualizationFunction) {
0244     return;
0245   }
0246   const MuonSegmentContainer& gotTruthSegs = m_inputTruthSegs(ctx);
0247   const std::string outputPath =
0248       std::format("HoughHistogram_{}_{}_{}_{}_{}.pdf",
0249                   MuonId::toString(bucketId.msStation()),
0250                   MuonId::toString(bucketId.side()), bucketId.sector(),
0251                   bucketId.measuresEta() ? "eta" : "phi", ctx.eventNumber);
0252   m_cfg.visualizationFunction(outputPath, bucketId, maxima, plane, axis,
0253                               gotTruthSegs, logger());
0254 }
0255 }  // namespace ActsExamples