TrackFinding/src/MuonHoughSeeder.cpp

0001 // This file is part of the Acts project.
0002 //
0003 // Copyright (C) 2020 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 http://mozilla.org/MPL/2.0/.
0008
0009 #include "ActsExamples/TrackFinding/MuonHoughSeeder.hpp"
0010
0011 #include "ActsExamples/EventData/MuonSimHit.hpp"
0012
0013 #include <algorithm>
0014 #include <cmath>
0015 #include <iterator>
0016 #include <ostream>
0017 #include <stdexcept>
0018 #include <variant>
0019
0020 #include "TBox.h"
0021 #include "TLatex.h"
0022 #include "TLegend.h"
0023
0024 namespace {
0025 /// map the station name integers in the CSV to the ATLAS station names
0026 static const std::map<int, std::string> stationDict{
0027     {0, "BIL"},  {1, "BIS"},  {7, "BIR"},  {2, "BML"},  {3, "BMS"},
0028     {8, "BMF"},  {53, "BME"}, {54, "BMG"}, {52, "BIM"}, {4, "BOL"},
0029     {5, "BOS"},  {9, "BOF"},  {10, "BOG"}, {6, "BEE"},  {14, "EEL"},
0030     {15, "EES"}, {13, "EIL"}, {49, "EIS"}, {17, "EML"}, {18, "EMS"},
0031     {20, "EOL"}, {21, "EOS"}};
0032
0033 }  // namespace
0034
0035 ActsExamples::MuonHoughSeeder::MuonHoughSeeder(
0036     ActsExamples::MuonHoughSeeder::Config cfg, Acts::Logging::Level lvl)
0037     : ActsExamples::IAlgorithm("MuonHoughSeeder", lvl),
0038       m_cfg(std::move(cfg)),
0039       m_logger(Acts::getDefaultLogger("MuonHoughSeeder", lvl)) {
0040   if (m_cfg.inDriftCircles.empty()) {
0041     throw std::invalid_argument(
0042         "MuonHoughSeeder: Missing drift circle collection");
0043   }
0044   if (m_cfg.inSimHits.empty()) {
0045     throw std::invalid_argument("MuonHoughSeeder: Missing sim hit collection");
0046   }
0047
0048   m_inputDriftCircles.initialize(m_cfg.inDriftCircles);
0049   m_inputSimHits.initialize(m_cfg.inSimHits);
0050 }
0051
0052 using PatternSeed = std::pair<double, double>;  // y0, tan theta
0053
0054 ActsExamples::ProcessCode ActsExamples::MuonHoughSeeder::execute(
0055     const AlgorithmContext& ctx) const {
0056   // read the hits and circles
0057   auto gotSH = m_inputSimHits(ctx);
0058   auto gotDC = m_inputDriftCircles(ctx);
0059
0060   // configure the binning of the hough plane
0061   Acts::HoughTransformUtils::HoughPlaneConfig planeCfg;
0062   planeCfg.nBinsX = 1000;
0063   planeCfg.nBinsY = 1000;
0064
0065   // instantiate the peak finder
0066   Acts::HoughTransformUtils::PeakFinders::IslandsAroundMaxConfig peakFinderCfg;
0067   peakFinderCfg.fractionCutoff = 0.7;
0068   peakFinderCfg.threshold = 3.;
0069   peakFinderCfg.minSpacingBetweenPeaks = {0., 30.};
0070
0071   // and map the hough plane to parameter ranges.
0072   // The first coordinate is tan(theta), the second is z0 in mm
0073   Acts::HoughTransformUtils::HoughAxisRanges axisRanges{-3., 3., -2000., 2000.};
0074
0075   // create the functions parametrising the hough space lines for drift circles.
0076   // Note that there are two solutions for each drift circle and angle
0077
0078   // left solution
0079   auto houghParam_fromDC_left = [](double tanTheta, const DriftCircle& DC) {
0080     return DC.y() - tanTheta * DC.z() -
0081            DC.rDrift() / std::cos(std::atan(tanTheta));
0082   };
0083   // right solution
0084   auto houghParam_fromDC_right = [](double tanTheta, const DriftCircle& DC) {
0085     return DC.y() - tanTheta * DC.z() +
0086            DC.rDrift() / std::cos(std::atan(tanTheta));
0087   };
0088
0089   // create the function parametrising the drift radius uncertainty
0090   auto houghWidth_fromDC = [](double, const DriftCircle& DC) {
0091     return std::min(DC.rDriftError() * 3.,
0092                     1.0);  // scale reported errors up to at least 1mm or 3
0093                            // times the reported error as drift circle calib not
0094                            // fully reliable at this stage
0095   };
0096
0097   // store the true parameters
0098   std::vector<PatternSeed> truePatterns;
0099
0100   // instantiate the hough plane
0101   Acts::HoughTransformUtils::HoughPlane<Acts::GeometryIdentifier::Value>
0102       houghPlane(planeCfg);
0103   // also insantiate the peak finder
0104   Acts::HoughTransformUtils::PeakFinders::IslandsAroundMax<
0105       Acts::GeometryIdentifier::Value>
0106       peakFinder(peakFinderCfg);
0107
0108   // loop pver true hirs
0109   for (auto& SH : gotSH) {
0110     // read the identifier
0111     muonMdtIdentifierFields detailedInfo =
0112         ActsExamples::splitId(SH.geometryId().value());
0113
0114     // store the true parameters
0115     truePatterns.emplace_back(SH.direction().y() / SH.direction().z(),
0116                               SH.fourPosition().y());
0117     // reset the hough plane
0118     houghPlane.reset();
0119     int foundDC = 0;
0120     // loop over drift circles
0121     for (DriftCircle& DC : gotDC) {
0122       if (DC.stationEta() == detailedInfo.stationEta &&
0123           DC.stationPhi() == detailedInfo.stationPhi &&
0124           DC.stationName() == detailedInfo.stationName) {
0125         // build a single identifier for the drift circles
0126         muonMdtIdentifierFields idf;
0127         idf.multilayer = DC.multilayer();
0128         idf.stationEta = DC.stationEta();
0129         idf.stationPhi = DC.stationPhi();
0130         idf.stationName = DC.stationName();
0131         idf.tubeLayer = DC.tubeLayer();
0132         idf.tube = DC.tube();
0133         auto identifier = compressId(idf);
0134         auto effectiveLayer = 3 * (DC.multilayer() - 1) + (DC.tubeLayer() - 1);
0135         ++foundDC;
0136         // populate the hough plane with both solutions.
0137         houghPlane.fill<DriftCircle>(DC, axisRanges, houghParam_fromDC_left,
0138                                      houghWidth_fromDC, identifier,
0139                                      effectiveLayer, 1.0);
0140         houghPlane.fill<DriftCircle>(DC, axisRanges, houghParam_fromDC_right,
0141                                      houghWidth_fromDC, identifier,
0142                                      effectiveLayer, 1.0);
0143       }
0144     }
0145     // now get the peaks
0146     auto maxima = peakFinder.findPeaks(houghPlane, axisRanges);
0147
0148     // visualisation in ROOT
0149     // represent the hough space as a TH2
0150     TH2D houghHistoForPlot("houghHist", "HoughPlane;tan(#theta);z0 [mm]",
0151                            houghPlane.nBinsX(), axisRanges.xMin,
0152                            axisRanges.xMax, houghPlane.nBinsY(),
0153                            axisRanges.yMin, axisRanges.yMax);
0154     for (int bx = 0; bx < houghHistoForPlot.GetNbinsX(); ++bx) {
0155       for (int by = 0; by < houghHistoForPlot.GetNbinsY(); ++by) {
0156         houghHistoForPlot.SetBinContent(bx + 1, by + 1,
0157                                         houghPlane.nHits(bx, by));
0158       }
0159     }
0160
0161     m_outCanvas->SetTitle(Form("Station %s, Eta %i, Phi %i",
0162                                stationDict.at(detailedInfo.stationName).c_str(),
0163                                (int)detailedInfo.stationEta,
0164                                (int)detailedInfo.stationPhi));
0165     houghHistoForPlot.SetTitle(
0166         Form("Station %s, Eta %i, Phi %i",
0167              stationDict.at(detailedInfo.stationName).c_str(),
0168              (int)detailedInfo.stationEta, (int)detailedInfo.stationPhi));
0169     m_outCanvas->cd();
0170     int maxHitsAsInt = static_cast<int>(houghPlane.maxHits());
0171     houghHistoForPlot.SetContour(maxHitsAsInt + 1);
0172     for (int k = 0; k < maxHitsAsInt + 1; ++k) {
0173       houghHistoForPlot.SetContourLevel(k, k - 0.5);
0174     }
0175     std::vector<std::unique_ptr<TMarker>> markers;
0176     std::vector<std::unique_ptr<TBox>> boxes;
0177     houghHistoForPlot.SetContourLevel(maxHitsAsInt + 1,
0178                                       houghHistoForPlot.GetMaximum() + 0.5);
0179     houghHistoForPlot.Draw("COLZ");
0180     // mark the true parameters
0181     auto trueMarker = std::make_unique<TMarker>(
0182         truePatterns.back().first, truePatterns.back().second, kOpenCrossX);
0183     trueMarker->SetMarkerSize(3);
0184     trueMarker->SetMarkerColor(kRed);
0185     trueMarker->Draw();
0186
0187     // now draw the hough maxima
0188     for (auto& max : maxima) {
0189       markers.push_back(std::make_unique<TMarker>(max.x, max.y, kFullSquare));
0190       markers.back()->SetMarkerSize(1);
0191       markers.back()->SetMarkerColor(kBlue);
0192       markers.back()->Draw();
0193       boxes.push_back(std::make_unique<TBox>(max.x - max.wx, max.y - max.wy,
0194                                              max.x + max.wx, max.y + max.wy));
0195       boxes.back()->SetLineColor(kBlue);
0196       boxes.back()->SetFillStyle(1001);
0197       boxes.back()->SetFillColorAlpha(kBlue, 0.1);
0198       boxes.back()->SetLineWidth(0);
0199       boxes.back()->Draw();
0200     }
0201     TLegend legend(0.5, 0.7, 1. - gPad->GetRightMargin(),
0202                    1. - gPad->GetTopMargin());
0203     legend.AddEntry(trueMarker.get(), "True coordinates");
0204     legend.SetBorderSize(0);
0205     legend.SetFillStyle(0);
0206     legend.Draw();
0207
0208     if (!boxes.empty()) {
0209       legend.AddEntry(markers.back().get(), "Hough maxima");
0210       legend.AddEntry(boxes.back().get(), "Hough uncertainties");
0211     }
0212     TLatex tl(gPad->GetLeftMargin() + 0.03, 1. - gPad->GetTopMargin() - 0.1,
0213               Form("%i drift circles on station", foundDC));
0214     tl.SetTextFont(43);
0215     tl.SetTextSize(24);
0216     tl.SetNDC();
0217     tl.Draw();
0218     m_outCanvas->SaveAs("HoughHistograms.pdf");
0219   }
0220
0221   ACTS_VERBOSE("SH: " << gotSH.size());
0222   ACTS_VERBOSE("DC: " << gotDC.size());
0223   return ActsExamples::ProcessCode::SUCCESS;
0224 }
0225
0226 ActsExamples::ProcessCode ActsExamples::MuonHoughSeeder::initialize() {
0227   // book the output canvas
0228   m_outCanvas = std::make_unique<TCanvas>("canvas", "", 800, 800);
0229   m_outCanvas->SaveAs("HoughHistograms.pdf[");
0230   m_outCanvas->SetRightMargin(0.12);
0231   m_outCanvas->SetLeftMargin(0.12);
0232   gStyle->SetPalette(kGreyScale);
0233   gStyle->SetOptStat(0);
0234   return ProcessCode::SUCCESS;
0235 }
0236 ActsExamples::ProcessCode ActsExamples::MuonHoughSeeder::finalize() {
0237   m_outCanvas->SaveAs("HoughHistograms.pdf]");
0238   return ProcessCode::SUCCESS;
0239 }