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     if (!m_cfg.extendWithPhi) {
0099       outMaxima.insert(outMaxima.end(),
0100                        std::make_move_iterator(etaMaxima.begin()),
0101                        std::make_move_iterator(etaMaxima.end()));
0102       continue;
0103     }
0104
0105     MuonHoughMaxContainer stationMaxima =
0106         extendMaximaWithPhi(ctx, std::move(etaMaxima), phiPlane);
0107     ACTS_VERBOSE(__func__ << "() - " << __LINE__ << " Extended the maxima to "
0108                           << stationMaxima.size() << " station maxima");
0109     outMaxima.insert(outMaxima.end(),
0110                      std::make_move_iterator(stationMaxima.begin()),
0111                      std::make_move_iterator(stationMaxima.end()));
0112   }
0113
0114   m_outputMaxima(ctx, std::move(outMaxima));
0115
0116   return ProcessCode::SUCCESS;
0117 }
0118
0119 MuonHoughMaxContainer MuonHoughSeeder::constructEtaMaxima(
0120     const AlgorithmContext& ctx, const MuonSpacePointBucket& bucket,
0121     HoughPlane_t& plane) const {
0122   MuonHoughMaxContainer etaMaxima{};
0123   AxisRange_t axisRanges{-3., 3., 100. * Acts::UnitConstants::m,
0124                          -100. * Acts::UnitConstants::m};
0125   /** brief Adapt the intercept axis ranges */
0126   for (const MuonSpacePoint& sp : bucket) {
0127     /** Require that the space point measures eta */
0128     if (!sp.id().measuresEta()) {
0129       continue;
0130     }
0131     const double y = sp.localPosition().y();
0132
0133     axisRanges.yMin = std::min(axisRanges.yMin, y - m_cfg.etaPlaneMarginIcept);
0134     axisRanges.yMax = std::max(axisRanges.yMax, y + m_cfg.etaPlaneMarginIcept);
0135   }
0136   /** Ranges are adapted. Now fill the hough plane */
0137   plane.reset();
0138   for (const MuonSpacePoint& sp : bucket) {
0139     /** Require that the space point measures eta */
0140     if (!sp.id().measuresEta()) {
0141       continue;
0142     }
0143     if (sp.id().technology() == MuonSpacePoint::MuonId::TechField::Mdt) {
0144       plane.fill<MuonSpacePoint>(sp, axisRanges, etaHoughParamDC_left,
0145                                  houghWidth_fromDC, &sp, sp.id().detLayer());
0146       plane.fill<MuonSpacePoint>(sp, axisRanges, etaHoughParamDC_right,
0147                                  houghWidth_fromDC, &sp, sp.id().detLayer());
0148     } else {
0149       plane.fill<MuonSpacePoint>(sp, axisRanges, etaHoughParam_strip,
0150                                  etaHoughWidth_strip, &sp, sp.id().detLayer());
0151     }
0152   }
0153   /** Extract the maxima from the peak */
0154   PeakFinderCfg_t peakFinderCfg{};
0155   peakFinderCfg.fractionCutoff = 0.7f;
0156   peakFinderCfg.threshold = 3.f;
0157   peakFinderCfg.minSpacingBetweenPeaks = {0.f, 30.f};
0158
0159   PeakFinder_t peakFinder{peakFinderCfg};
0160
0161   /// Fetch the first set of peaks
0162   const MaximumVec_t maxima = peakFinder.findPeaks(plane, axisRanges);
0163   if (maxima.empty()) {
0164     ACTS_VERBOSE(__func__ << "() - " << __LINE__
0165                           << " No maxima found for bucket");
0166     return etaMaxima;
0167   }
0168   /** Create a measurement vector of the pure phi hits to be copied on all eta
0169    * maxima */
0170
0171   MuonHoughMaximum::HitVec phiHits{};
0172   if (m_cfg.extendWithPhi) {
0173     for (const MuonSpacePoint& sp : bucket) {
0174       if (!sp.id().measuresEta()) {
0175         phiHits.emplace_back(&sp);
0176       }
0177     }
0178   }
0179   for (const Maximum_t& peakMax : maxima) {
0180     MuonHoughMaximum::HitVec hits{peakMax.hitIdentifiers.begin(),
0181                                   peakMax.hitIdentifiers.end()};
0182     hits.insert(hits.end(), phiHits.begin(), phiHits.end());
0183     etaMaxima.emplace_back(peakMax.x, peakMax.y, std::move(hits));
0184   }
0185   MuonId bucketId{bucket.front().id()};
0186   bucketId.setCoordFlags(true, false);
0187   displayMaxima(ctx, bucketId, maxima, plane, axisRanges);
0188   return etaMaxima;
0189 }
0190 MuonHoughMaxContainer MuonHoughSeeder::extendMaximaWithPhi(
0191     const AlgorithmContext& ctx, MuonHoughMaxContainer&& etaMaxima,
0192     HoughPlane_t& plane) const {
0193   MuonHoughMaxContainer outMaxima{};
0194
0195   PeakFinderCfg_t peakFinderCfg{};
0196   peakFinderCfg.fractionCutoff = 0.7;
0197   peakFinderCfg.threshold = 2.;
0198   peakFinderCfg.minSpacingBetweenPeaks = {0., 30.};
0199
0200   PeakFinder_t peakFinder{peakFinderCfg};
0201
0202   for (MuonHoughMaximum& etaMax : etaMaxima) {
0203     plane.reset();
0204     AxisRange_t axisRanges{-3., 3., 100. * Acts::UnitConstants::m,
0205                            -100. * Acts::UnitConstants::m};
0206
0207     unsigned nPhiHits{0};
0208     for (const MuonSpacePoint* sp : etaMax.hits()) {
0209       if (!sp->id().measuresPhi()) {
0210         continue;
0211       }
0212       ++nPhiHits;
0213       const double x = sp->localPosition().x();
0214       axisRanges.yMin =
0215           std::min(axisRanges.yMin, x - m_cfg.phiPlaneMarginIcept);
0216       axisRanges.yMax =
0217           std::max(axisRanges.yMax, x + m_cfg.phiPlaneMarginIcept);
0218     }
0219     if (nPhiHits < 2) {
0220       outMaxima.emplace_back(std::move(etaMax));
0221       continue;
0222     }
0223     MuonHoughMaximum::HitVec etaHits{};
0224     etaHits.reserve(etaMax.hits().size());
0225     for (const MuonSpacePoint* sp : etaMax.hits()) {
0226       if (!sp->id().measuresPhi()) {
0227         etaHits.push_back(sp);
0228         continue;
0229       }
0230       plane.fill<MuonSpacePoint>(*sp, axisRanges, phiHoughParam_strip,
0231                                  phiHoughWidth_strip, sp, sp->id().detLayer());
0232     }
0233     const MaximumVec_t phiMaxima = peakFinder.findPeaks(plane, axisRanges);
0234     for (const Maximum_t& max : phiMaxima) {
0235       MuonHoughMaximum::HitVec hits{max.hitIdentifiers.begin(),
0236                                     max.hitIdentifiers.end()};
0237       hits.insert(hits.end(), etaHits.begin(), etaHits.end());
0238       outMaxima.emplace_back(max.x, max.y, etaMax.tanBeta(),
0239                              etaMax.interceptY(), std::move(hits));
0240     }
0241     MuonId bucketId{etaMax.hits().front()->id()};
0242     bucketId.setCoordFlags(false, true);
0243     displayMaxima(ctx, bucketId, phiMaxima, plane, axisRanges);
0244     if (phiMaxima.empty()) {
0245       outMaxima.emplace_back(std::move(etaMax));
0246     }
0247   }
0248   return outMaxima;
0249 }
0250
0251 void MuonHoughSeeder::displayMaxima(const AlgorithmContext& ctx,
0252                                     const MuonId& bucketId,
0253                                     const MaximumVec_t& maxima,
0254                                     const HoughPlane_t& plane,
0255                                     const AxisRange_t& axis) const {
0256   if (!m_cfg.dumpVisualization || !m_cfg.visualizationFunction) {
0257     return;
0258   }
0259   const MuonSegmentContainer& gotTruthSegs = m_inputTruthSegs(ctx);
0260   const std::string outputPath =
0261       std::format("HoughHistogram_{}_{}_{}_{}_{}.pdf",
0262                   MuonId::toString(bucketId.msStation()),
0263                   MuonId::toString(bucketId.side()), bucketId.sector(),
0264                   bucketId.measuresEta() ? "eta" : "phi", ctx.eventNumber);
0265   m_cfg.visualizationFunction(outputPath, bucketId, maxima, plane, axis,
0266                               gotTruthSegs, logger());
0267 }
0268 }  // namespace ActsExamples