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 "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     // scale reported errors up to at least 1mm or 3 times the reported error as
0092     // drift circle calib not fully reliable at this stage
0093     return std::min(DC.rDriftError() * 3., 1.0);
0094   };
0095
0096   // store the true parameters
0097   std::vector<PatternSeed> truePatterns;
0098
0099   // instantiate the hough plane
0100   Acts::HoughTransformUtils::HoughPlane<Acts::GeometryIdentifier::Value>
0101       houghPlane(planeCfg);
0102   // also instantiate the peak finder
0103   Acts::HoughTransformUtils::PeakFinders::IslandsAroundMax<
0104       Acts::GeometryIdentifier::Value>
0105       peakFinder(peakFinderCfg);
0106
0107   // loop over true hits
0108   for (auto& SH : gotSH) {
0109     // read the identifier
0110     MuonMdtIdentifierFields detailedInfo =
0111         ActsExamples::splitId(SH.geometryId().value());
0112     // store the true parameters
0113     truePatterns.emplace_back(SH.direction().y() / SH.direction().z(),
0114                               SH.fourPosition().y());
0115     // ACTS_VERBOSE("station name=" << static_cast<int>(SH.stationName));
0116     ACTS_VERBOSE("direction = " << SH.direction().y());
0117     ACTS_VERBOSE("fourposition y = " << SH.fourPosition().y());
0118     std::cin.ignore();
0119     // reset the hough plane
0120     houghPlane.reset();
0121     int foundDC = 0;
0122     // loop over drift circles
0123     for (DriftCircle& DC : gotDC) {
0124       if (DC.stationEta() == detailedInfo.stationEta &&
0125           DC.stationPhi() == detailedInfo.stationPhi &&
0126           DC.stationName() == detailedInfo.stationName) {
0127         // build a single identifier for the drift circles
0128         MuonMdtIdentifierFields idf;
0129         idf.multilayer = DC.multilayer();
0130         idf.stationEta = DC.stationEta();
0131         idf.stationPhi = DC.stationPhi();
0132         idf.stationName = DC.stationName();
0133         idf.tubeLayer = DC.tubeLayer();
0134         idf.tube = DC.tube();
0135         auto identifier = compressId(idf);
0136         auto effectiveLayer = 3 * (DC.multilayer() - 1) + (DC.tubeLayer() - 1);
0137         ++foundDC;
0138         // populate the hough plane with both solutions.
0139         houghPlane.fill<DriftCircle>(DC, axisRanges, houghParam_fromDC_left,
0140                                      houghWidth_fromDC, identifier,
0141                                      effectiveLayer, 1.0);
0142         houghPlane.fill<DriftCircle>(DC, axisRanges, houghParam_fromDC_right,
0143                                      houghWidth_fromDC, identifier,
0144                                      effectiveLayer, 1.0);
0145       }
0146     }
0147     // now get the peaks
0148     auto maxima = peakFinder.findPeaks(houghPlane, axisRanges);
0149
0150     // visualisation in ROOT
0151     // represent the hough space as a TH2
0152     TH2D houghHistoForPlot("houghHist", "HoughPlane;tan(#theta);z0 [mm]",
0153                            houghPlane.nBinsX(), axisRanges.xMin,
0154                            axisRanges.xMax, houghPlane.nBinsY(),
0155                            axisRanges.yMin, axisRanges.yMax);
0156     for (int bx = 0; bx < houghHistoForPlot.GetNbinsX(); ++bx) {
0157       for (int by = 0; by < houghHistoForPlot.GetNbinsY(); ++by) {
0158         houghHistoForPlot.SetBinContent(bx + 1, by + 1,
0159                                         houghPlane.nHits(bx, by));
0160       }
0161     }
0162     m_outCanvas->SetTitle(Form("Station %s, Eta %i, Phi %i",
0163                                stationDict.at(detailedInfo.stationName).c_str(),
0164                                static_cast<int>(detailedInfo.stationEta),
0165                                static_cast<int>(detailedInfo.stationPhi)));
0166     houghHistoForPlot.SetTitle(
0167         Form("Station %s, Eta %i, Phi %i",
0168              stationDict.at(detailedInfo.stationName).c_str(),
0169              static_cast<int>(detailedInfo.stationEta),
0170              static_cast<int>(detailedInfo.stationPhi)));
0171     m_outCanvas->cd();
0172     int maxHitsAsInt = static_cast<int>(houghPlane.maxHits());
0173     houghHistoForPlot.SetContour(maxHitsAsInt + 1);
0174     for (int k = 0; k < maxHitsAsInt + 1; ++k) {
0175       houghHistoForPlot.SetContourLevel(k, k - 0.5);
0176     }
0177     std::vector<std::unique_ptr<TMarker>> markers;
0178     std::vector<std::unique_ptr<TBox>> boxes;
0179     houghHistoForPlot.SetContourLevel(maxHitsAsInt + 1,
0180                                       houghHistoForPlot.GetMaximum() + 0.5);
0181     houghHistoForPlot.Draw("COLZ");
0182     // mark the true parameters
0183     auto trueMarker = std::make_unique<TMarker>(
0184         truePatterns.back().first, truePatterns.back().second, kOpenCrossX);
0185     trueMarker->SetMarkerSize(3);
0186     trueMarker->SetMarkerColor(kRed);
0187     trueMarker->Draw();
0188
0189     // now draw the hough maxima
0190     for (auto& max : maxima) {
0191       markers.push_back(std::make_unique<TMarker>(max.x, max.y, kFullSquare));
0192       markers.back()->SetMarkerSize(1);
0193       markers.back()->SetMarkerColor(kBlue);
0194       markers.back()->Draw();
0195       boxes.push_back(std::make_unique<TBox>(max.x - max.wx, max.y - max.wy,
0196                                              max.x + max.wx, max.y + max.wy));
0197       boxes.back()->SetLineColor(kBlue);
0198       boxes.back()->SetFillStyle(1001);
0199       boxes.back()->SetFillColorAlpha(kBlue, 0.1);
0200       boxes.back()->SetLineWidth(0);
0201       boxes.back()->Draw();
0202     }
0203     TLegend legend(0.5, 0.7, 1. - gPad->GetRightMargin(),
0204                    1. - gPad->GetTopMargin());
0205     legend.AddEntry(trueMarker.get(), "True coordinates");
0206     legend.SetBorderSize(0);
0207     legend.SetFillStyle(0);
0208     legend.Draw();
0209
0210     if (!boxes.empty()) {
0211       legend.AddEntry(markers.back().get(), "Hough maxima");
0212       legend.AddEntry(boxes.back().get(), "Hough uncertainties");
0213     }
0214     TLatex tl(gPad->GetLeftMargin() + 0.03, 1. - gPad->GetTopMargin() - 0.1,
0215               Form("%i drift circles on station", foundDC));
0216     tl.SetTextFont(43);
0217     tl.SetTextSize(24);
0218     tl.SetNDC();
0219     tl.Draw();
0220     m_outCanvas->SaveAs("HoughHistograms.pdf");
0221   }
0222
0223   ACTS_VERBOSE("SH: " << gotSH.size());
0224   ACTS_VERBOSE("DC: " << gotDC.size());
0225   return ActsExamples::ProcessCode::SUCCESS;
0226 }
0227
0228 ActsExamples::ProcessCode ActsExamples::MuonHoughSeeder::initialize() {
0229   // book the output canvas
0230   m_outCanvas = std::make_unique<TCanvas>("canvas", "", 800, 800);
0231   m_outCanvas->SaveAs("HoughHistograms.pdf[");
0232   m_outCanvas->SetRightMargin(0.12);
0233   m_outCanvas->SetLeftMargin(0.12);
0234   gStyle->SetPalette(kGreyScale);
0235   gStyle->SetOptStat(0);
0236   return ProcessCode::SUCCESS;
0237 }
0238 ActsExamples::ProcessCode ActsExamples::MuonHoughSeeder::finalize() {
0239   m_outCanvas->SaveAs("HoughHistograms.pdf]");
0240   return ProcessCode::SUCCESS;
0241 }