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