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 #include "TBox.h"
0022 #include "TCanvas.h"
0023 #include "TH2D.h"
0024 #include "TLegend.h"
0025 #include "TMarker.h"
0026 #include "TStyle.h"
0027
0028 namespace ActsExamples {
0029
0030 // left solution
0031 auto etaHoughParamDC_left = [](double tanAlpha, const MuonSpacePoint& DC) {
0032   return DC.localPosition().y() - tanAlpha * DC.localPosition().z() -
0033          DC.driftRadius() * std::sqrt(1. + tanAlpha * tanAlpha);
0034 };
0035 // right solution
0036 auto etaHoughParamDC_right = [](double tanAlpha, const MuonSpacePoint& DC) {
0037   return DC.localPosition().y() - tanAlpha * DC.localPosition().z() +
0038          DC.driftRadius() * std::sqrt(1. + tanAlpha * tanAlpha);
0039 };
0040
0041 // create the function parametrising the drift radius uncertainty
0042 auto houghWidth_fromDC = [](double, const MuonSpacePoint& DC) {
0043   // scale reported errors up to at least 1mm or 3 times the reported error as
0044   // drift circle calib not fully reliable at this stage
0045   return std::min(std::sqrt(DC.covariance()[Acts::eY]) * 3., 1.0);
0046 };
0047
0048 /// strip solution
0049 auto etaHoughParam_strip = [](double tanAlpha, const MuonSpacePoint& strip) {
0050   return strip.localPosition().y() - tanAlpha * strip.localPosition().z();
0051 };
0052 auto phiHoughParam_strip = [](double tanBeta, const MuonSpacePoint& strip) {
0053   return strip.localPosition().x() - tanBeta * strip.localPosition().z();
0054 };
0055
0056 /// @brief Strip uncertainty
0057 auto etaHoughWidth_strip = [](double, const MuonSpacePoint& strip) {
0058   return std::sqrt(strip.covariance()[Acts::eY]) * 3.;
0059 };
0060 auto phiHoughWidth_strip = [](double, const MuonSpacePoint& strip) {
0061   return std::sqrt(strip.covariance()[Acts::eX]) * 3.;
0062 };
0063
0064 MuonHoughSeeder::MuonHoughSeeder(MuonHoughSeeder::Config cfg,
0065                                  Acts::Logging::Level lvl)
0066     : IAlgorithm("MuonHoughSeeder", lvl),
0067       m_cfg(std::move(cfg)),
0068       m_logger(Acts::getDefaultLogger(name(), lvl)) {
0069   if (m_cfg.inSpacePoints.empty()) {
0070     throw std::invalid_argument(
0071         "MuonHoughSeeder: Missing drift circle collection");
0072   }
0073   if (m_cfg.inTruthSegments.empty()) {
0074     throw std::invalid_argument(
0075         "MuonHoughSeeder: Missing truth segment collection");
0076   }
0077
0078   m_inputSpacePoints.initialize(m_cfg.inSpacePoints);
0079   m_inputTruthSegs.initialize(m_cfg.inTruthSegments);
0080   m_outputMaxima.initialize(m_cfg.outHoughMax);
0081 }
0082
0083 MuonHoughSeeder::~MuonHoughSeeder() = default;
0084
0085 ProcessCode MuonHoughSeeder::execute(const AlgorithmContext& ctx) const {
0086   // read the hits and circles
0087   const MuonSpacePointContainer& gotSpacePoints = m_inputSpacePoints(ctx);
0088
0089   // configure the binning of the hough plane
0090   Acts::HoughTransformUtils::HoughPlaneConfig etaPlaneCfg;
0091   etaPlaneCfg.nBinsX = m_cfg.nBinsTanTheta;
0092   etaPlaneCfg.nBinsY = m_cfg.nBinsY0;
0093
0094   Acts::HoughTransformUtils::HoughPlaneConfig phiPlaneCfg;
0095   phiPlaneCfg.nBinsX = m_cfg.nBinsTanPhi;
0096   phiPlaneCfg.nBinsY = m_cfg.nBinsX0;
0097
0098   // instantiate the hough plane
0099   HoughPlane_t etaPlane{etaPlaneCfg};
0100   HoughPlane_t phiPlane{phiPlaneCfg};
0101   MuonHoughMaxContainer outMaxima{};
0102
0103   for (const MuonSpacePointContainer::value_type& bucket : gotSpacePoints) {
0104     MuonHoughMaxContainer etaMaxima = constructEtaMaxima(ctx, bucket, etaPlane);
0105     ACTS_VERBOSE(__func__ << "() - " << __LINE__ << " Found "
0106                           << etaMaxima.size() << " eta maxima");
0107     MuonHoughMaxContainer stationMaxima =
0108         extendMaximaWithPhi(ctx, std::move(etaMaxima), phiPlane);
0109     ACTS_VERBOSE(__func__ << "() - " << __LINE__ << " Extended the maxima to "
0110                           << stationMaxima.size() << " station maxima");
0111     outMaxima.insert(outMaxima.end(),
0112                      std::make_move_iterator(stationMaxima.begin()),
0113                      std::make_move_iterator(stationMaxima.end()));
0114   }
0115
0116   m_outputMaxima(ctx, std::move(outMaxima));
0117
0118   return ProcessCode::SUCCESS;
0119 }
0120
0121 MuonHoughMaxContainer MuonHoughSeeder::constructEtaMaxima(
0122     const AlgorithmContext& ctx, const MuonSpacePointBucket& bucket,
0123     HoughPlane_t& plane) const {
0124   MuonHoughMaxContainer etaMaxima{};
0125   AxisRange_t axisRanges{-3., 3., 100. * Acts::UnitConstants::m,
0126                          -100. * Acts::UnitConstants::m};
0127   /** brief Adapt the intercept axis ranges */
0128   for (const MuonSpacePoint& sp : bucket) {
0129     /** Require that the space point measures eta */
0130     if (!sp.id().measuresEta()) {
0131       continue;
0132     }
0133     const double y = sp.localPosition().y();
0134
0135     axisRanges.yMin = std::min(axisRanges.yMin, y - m_cfg.etaPlaneMarginIcept);
0136     axisRanges.yMax = std::min(axisRanges.yMax, y + m_cfg.etaPlaneMarginIcept);
0137   }
0138   /** Ranges are adapted. Now fill the hough plane */
0139   plane.reset();
0140   for (const MuonSpacePoint& sp : bucket) {
0141     /** Require that the space point measures eta */
0142     if (!sp.id().measuresEta()) {
0143       continue;
0144     }
0145     if (sp.id().technology() == MuonSpacePoint::MuonId::TechField::Mdt) {
0146       plane.fill<MuonSpacePoint>(sp, axisRanges, etaHoughParamDC_left,
0147                                  houghWidth_fromDC, &sp, sp.id().detLayer());
0148       plane.fill<MuonSpacePoint>(sp, axisRanges, etaHoughParamDC_right,
0149                                  houghWidth_fromDC, &sp, sp.id().detLayer());
0150     } else {
0151       plane.fill<MuonSpacePoint>(sp, axisRanges, etaHoughParam_strip,
0152                                  etaHoughWidth_strip, &sp, sp.id().detLayer());
0153     }
0154   }
0155   /** Extract the maxima from the peak */
0156   PeakFinderCfg_t peakFinderCfg{};
0157   peakFinderCfg.fractionCutoff = 0.7f;
0158   peakFinderCfg.threshold = 3.f;
0159   peakFinderCfg.minSpacingBetweenPeaks = {0.f, 30.f};
0160
0161   PeakFinder_t peakFinder{peakFinderCfg};
0162
0163   /// Fetch the first set of peaks
0164   const MaximumVec_t maxima = peakFinder.findPeaks(plane, axisRanges);
0165   if (maxima.empty()) {
0166     ACTS_VERBOSE(__func__ << "() - " << __LINE__
0167                           << " No maxima found for bucket");
0168     return etaMaxima;
0169   }
0170   /** Create a measurement vector of the pure phi hits to be copied on all eta
0171    * maxima */
0172   MuonHoughMaximum::HitVec phiHits{};
0173   for (const MuonSpacePoint& sp : bucket) {
0174     if (!sp.id().measuresEta()) {
0175       phiHits.emplace_back(&sp);
0176     }
0177   }
0178   for (const Maximum_t& peakMax : maxima) {
0179     MuonHoughMaximum::HitVec hits{peakMax.hitIdentifiers.begin(),
0180                                   peakMax.hitIdentifiers.end()};
0181     hits.insert(hits.end(), phiHits.begin(), phiHits.end());
0182     etaMaxima.emplace_back(peakMax.x, peakMax.y, std::move(hits));
0183   }
0184   MuonId bucketId{bucket.front().id()};
0185   bucketId.setCoordFlags(true, false);
0186   displayMaxima(ctx, bucketId, maxima, plane, axisRanges);
0187   return etaMaxima;
0188 }
0189 MuonHoughMaxContainer MuonHoughSeeder::extendMaximaWithPhi(
0190     const AlgorithmContext& ctx, MuonHoughMaxContainer&& etaMaxima,
0191     HoughPlane_t& plane) const {
0192   MuonHoughMaxContainer outMaxima{};
0193
0194   PeakFinderCfg_t peakFinderCfg{};
0195   peakFinderCfg.fractionCutoff = 0.7;
0196   peakFinderCfg.threshold = 2.;
0197   peakFinderCfg.minSpacingBetweenPeaks = {0., 30.};
0198
0199   PeakFinder_t peakFinder{peakFinderCfg};
0200
0201   for (MuonHoughMaximum& etaMax : etaMaxima) {
0202     plane.reset();
0203     AxisRange_t axisRanges{-3., 3., 100. * Acts::UnitConstants::m,
0204                            -100. * Acts::UnitConstants::m};
0205
0206     unsigned nPhiHits{0};
0207     for (const MuonSpacePoint* sp : etaMax.hits()) {
0208       if (!sp->id().measuresPhi()) {
0209         continue;
0210       }
0211       ++nPhiHits;
0212       const double x = sp->localPosition().x();
0213       axisRanges.yMax =
0214           std::min(axisRanges.yMax, x + m_cfg.phiPlaneMarginIcept);
0215       axisRanges.yMin =
0216           std::min(axisRanges.yMin, x - m_cfg.phiPlaneMarginIcept);
0217     }
0218     if (nPhiHits < 2) {
0219       outMaxima.emplace_back(std::move(etaMax));
0220       continue;
0221     }
0222     MuonHoughMaximum::HitVec etaHits{};
0223     etaHits.reserve(etaMax.hits().size());
0224     for (const MuonSpacePoint* sp : etaMax.hits()) {
0225       if (!sp->id().measuresPhi()) {
0226         etaHits.push_back(sp);
0227         continue;
0228       }
0229       plane.fill<MuonSpacePoint>(*sp, axisRanges, phiHoughParam_strip,
0230                                  phiHoughWidth_strip, sp, sp->id().detLayer());
0231     }
0232     const MaximumVec_t phiMaxima = peakFinder.findPeaks(plane, axisRanges);
0233     for (const Maximum_t& max : phiMaxima) {
0234       MuonHoughMaximum::HitVec hits{max.hitIdentifiers.begin(),
0235                                     max.hitIdentifiers.end()};
0236       hits.insert(hits.end(), etaHits.begin(), etaHits.end());
0237       outMaxima.emplace_back(max.x, max.y, etaMax.tanBeta(),
0238                              etaMax.interceptY(), std::move(hits));
0239     }
0240     MuonId bucketId{etaMax.hits().front()->id()};
0241     bucketId.setCoordFlags(false, true);
0242     displayMaxima(ctx, bucketId, phiMaxima, plane, axisRanges);
0243   }
0244   return outMaxima;
0245 }
0246
0247 void MuonHoughSeeder::displayMaxima(const AlgorithmContext& ctx,
0248                                     const MuonId& bucketId,
0249                                     const MaximumVec_t& maxima,
0250                                     const HoughPlane_t& plane,
0251                                     const AxisRange_t& axis) const {
0252   if (!m_outCanvas) {
0253     return;
0254   }
0255   static std::mutex canvasMutex{};
0256   // read the hits and circles
0257   const MuonSegmentContainer& gotTruthSegs = m_inputTruthSegs(ctx);
0258   std::lock_guard guard{canvasMutex};
0259   std::vector<std::unique_ptr<TObject>> primitives;
0260
0261   /// Save the hough accumulator as histogram
0262   TH2D houghHistoForPlot("houghHist", "HoughPlane;tan(#alpha);z0 [mm]",
0263                          plane.nBinsX(), axis.xMin, axis.xMax, plane.nBinsY(),
0264                          axis.yMin, axis.yMax);
0265   houghHistoForPlot.SetTitle(std::format("Station {:}, side {:}, sector {:2d}",
0266                                          MuonId::toString(bucketId.msStation()),
0267                                          MuonId::toString(bucketId.side()),
0268                                          bucketId.sector())
0269                                  .c_str());
0270
0271   /** Copy the plane content into the histogram */
0272   for (int bx = 0; bx < houghHistoForPlot.GetNbinsX(); ++bx) {
0273     for (int by = 0; by < houghHistoForPlot.GetNbinsY(); ++by) {
0274       houghHistoForPlot.SetBinContent(bx + 1, by + 1, plane.nHits(bx, by));
0275     }
0276   }
0277   /** Set the contours */
0278   auto maxHitsAsInt = static_cast<int>(plane.maxHits());
0279   houghHistoForPlot.SetContour(maxHitsAsInt + 1);
0280   for (int k = 0; k < maxHitsAsInt + 1; ++k) {
0281     houghHistoForPlot.SetContourLevel(k, k - 0.5);
0282   }
0283
0284   auto legend = std::make_unique<TLegend>(0.5, 0.7, 1. - gPad->GetRightMargin(),
0285                                           1. - gPad->GetTopMargin());
0286   legend->SetBorderSize(0);
0287   legend->SetFillStyle(0);
0288
0289   /** Fetch the true parameters */
0290   MuonSegmentContainer::const_iterator truthItr = gotTruthSegs.begin();
0291   const int trueCoord = bucketId.measuresEta() ? Acts::eY : Acts::eX;
0292   while ((truthItr = std::find_if(truthItr, gotTruthSegs.end(),
0293                                   [bucketId](const MuonSegment& seg) {
0294                                     return seg.id().sameStation(bucketId);
0295                                   })) != gotTruthSegs.end()) {
0296     const MuonSegment& truthSeg{*truthItr};
0297     const float tanAlpha =
0298         truthSeg.localDirection()[trueCoord] / truthSeg.localDirection().z();
0299     const float intercept = truthSeg.localPosition()[trueCoord];
0300
0301     auto trueMarker =
0302         std::make_unique<TMarker>(tanAlpha, intercept, kOpenCrossX);
0303     trueMarker->SetMarkerSize(3);
0304     trueMarker->SetMarkerColor(kRed);
0305     legend->AddEntry(trueMarker.get(), "True coordinates");
0306     primitives.push_back(std::move(trueMarker));
0307     ACTS_VERBOSE("Draw true segment " << truthSeg.id()
0308                                       << " in bucket: " << bucketId);
0309     ++truthItr;
0310   }
0311
0312   bool addedLeg{false};
0313   for (const Maximum_t& max : maxima) {
0314     auto marker = std::make_unique<TMarker>(max.x, max.y, kFullSquare);
0315     marker->SetMarkerSize(1);
0316     marker->SetMarkerColor(kBlue);
0317
0318     auto box = std::make_unique<TBox>(max.x - max.wx, max.y - max.wy,
0319                                       max.x + max.wx, max.y + max.wy);
0320     box->SetLineColor(kBlue);
0321     box->SetFillStyle(1001);
0322     box->SetFillColorAlpha(kBlue, 0.1);
0323     box->SetLineWidth(0);
0324     if (!addedLeg) {
0325       legend->AddEntry(marker.get(), "Hough maxima");
0326       legend->AddEntry(box.get(), "Hough uncertainties");
0327       addedLeg = true;
0328     }
0329     primitives.push_back(std::move(box));
0330     primitives.push_back(std::move(marker));
0331   }
0332   primitives.emplace_back(std::move(legend));
0333   for (auto& prim : primitives) {
0334     prim->Draw();
0335   }
0336
0337   m_outCanvas->SaveAs("HoughHistograms.pdf");
0338 }
0339 ProcessCode MuonHoughSeeder::initialize() {
0340   // book the output canvas
0341   m_outCanvas = std::make_unique<TCanvas>("canvas", "", 800, 800);
0342   m_outCanvas->SaveAs("HoughHistograms.pdf[");
0343   m_outCanvas->SetRightMargin(0.12);
0344   m_outCanvas->SetLeftMargin(0.12);
0345   gStyle->SetPalette(kGreyScale);
0346   gStyle->SetOptStat(0);
0347   return ProcessCode::SUCCESS;
0348 }
0349 ProcessCode MuonHoughSeeder::finalize() {
0350   m_outCanvas->SaveAs("HoughHistograms.pdf]");
0351   return ProcessCode::SUCCESS;
0352 }
0353 }  // namespace ActsExamples