Algorithms/TrackFinding/MuonHoughSeederTests.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 <boost/test/unit_test.hpp>
0010
0011 #include "Acts/Definitions/Algebra.hpp"
0012 #include "Acts/Definitions/Units.hpp"
0013 #include "Acts/Geometry/GeometryIdentifier.hpp"
0014 #include "Acts/Utilities/Logger.hpp"
0015 #include "ActsExamples/EventData/MuonHoughMaximum.hpp"
0016 #include "ActsExamples/EventData/MuonSegment.hpp"
0017 #include "ActsExamples/EventData/MuonSpacePoint.hpp"
0018 #include "ActsExamples/Framework/AlgorithmContext.hpp"
0019 #include "ActsExamples/Framework/DataHandle.hpp"
0020 #include "ActsExamples/Framework/WhiteBoard.hpp"
0021 #include "ActsExamples/TrackFinding/MuonHoughSeeder.hpp"
0022
0023 #include <array>
0024 #include <cmath>
0025 #include <cstddef>
0026 #include <cstdint>
0027 #include <limits>
0028 #include <vector>
0029
0030 namespace ActsTests {
0031
0032 namespace {
0033
0034 struct DriftCircle {
0035   double y = 0.0;
0036   double z = 0.0;
0037   double rDrift = 0.0;
0038   double rDriftError = 0.0;
0039 };
0040
0041 /// @brief Prepare MuonId
0042 ActsExamples::MuonSpacePoint::MuonId makeMdtEtaId(std::uint8_t layer,
0043                                                   std::uint16_t channel) {
0044   using MuonId = ActsExamples::MuonSpacePoint::MuonId;
0045
0046   MuonId id{};
0047   id.setChamber(MuonId::StationName::BIL, MuonId::DetSide::A, 1u,
0048                 MuonId::TechField::Mdt);
0049   id.setLayAndCh(layer, channel);
0050   id.setCoordFlags(true, false, false);
0051
0052   return id;
0053 }
0054
0055 /// @brief Prepare example MuonSpaceContainer with one bucket based on test case from HoughTransformUtilsTests.cpp
0056 ActsExamples::MuonSpacePointContainer makeDriftCircleSpacePoints() {
0057   constexpr double uncert = 0.3;
0058
0059   const std::array<DriftCircle, 6> driftCircles{
0060       DriftCircle{-427.981, -225.541, 14.5202, uncert},
0061       DriftCircle{-412.964, -199.530, 1.66237, uncert},
0062       DriftCircle{-427.981, -173.519, 12.3176, uncert},
0063       DriftCircle{-427.981, 173.519, 1.5412, uncert},
0064       DriftCircle{-442.999, 199.530, 12.3937, uncert},
0065       DriftCircle{-427.981, 225.541, 3.77967, uncert},
0066   };
0067
0068   ActsExamples::MuonSpacePointContainer spacePoints{};
0069   spacePoints.emplace_back();
0070
0071   auto& bucket = spacePoints.back();
0072   bucket.reserve(driftCircles.size());
0073
0074   for (std::size_t i = 0; i < driftCircles.size(); ++i) {
0075     const DriftCircle& dc = driftCircles[i];
0076
0077     ActsExamples::MuonSpacePoint& sp = bucket.emplace_back();
0078
0079     sp.setGeometryId(Acts::GeometryIdentifier{i + 1u});
0080     sp.setId(makeMdtEtaId(static_cast<std::uint8_t>(i + 1u),
0081                           static_cast<std::uint16_t>(i + 1u)));
0082
0083     sp.defineCoordinates(Acts::Vector3{0.0, dc.y, dc.z}, Acts::Vector3::UnitX(),
0084                          Acts::Vector3::UnitY());
0085
0086     sp.setRadius(dc.rDrift);
0087     sp.setTime(0.0);
0088     sp.setCovariance(0.0, dc.rDriftError * dc.rDriftError, 0.0);
0089   }
0090
0091   return spacePoints;
0092 }
0093
0094 }  // namespace
0095
0096 BOOST_AUTO_TEST_SUITE(MuonHoughTransformSuite)
0097
0098 BOOST_AUTO_TEST_CASE(muon_hough_seeder_drift_circle_sanity) {
0099   // Truth from the original HoughTransformUtils unit test.
0100   constexpr double trueTanTheta = -0.0401472 / 0.994974;
0101   constexpr double trueInterceptY = -422.612;
0102
0103   ActsExamples::MuonHoughSeeder::Config cfg{};
0104   cfg.inSpacePoints = "MuonSpacePoints";
0105   cfg.inTruthSegments =
0106       "TruthSegments";  // this is required even if empty in ctx
0107   cfg.outHoughMax = "MuonHoughMaxima";
0108
0109   cfg.nBinsTanTheta = 1000;
0110   cfg.nBinsY0 = 1000;
0111
0112   cfg.nBinsTanPhi = 10;
0113   cfg.nBinsX0 = 10;
0114
0115   cfg.etaPlaneMarginIcept = 2.0 * Acts::UnitConstants::m;
0116   cfg.phiPlaneMarginIcept = 2.0 * Acts::UnitConstants::m;
0117
0118   cfg.dumpVisualization = false;
0119
0120   ActsExamples::MuonHoughSeeder seeder{
0121       cfg, Acts::getDefaultLogger("MuonHoughSeederTest", Acts::Logging::INFO)};
0122
0123   ActsExamples::WhiteBoard eventStore{};
0124   ActsExamples::AlgorithmContext ctx{0, 0, eventStore, 0};
0125
0126   ActsExamples::WriteDataHandle<ActsExamples::MuonSpacePointContainer>
0127       spacePointHandle{&seeder, "TestInputSpacePoints"};
0128   spacePointHandle.initialize(cfg.inSpacePoints);
0129   spacePointHandle(ctx, makeDriftCircleSpacePoints());
0130
0131   BOOST_REQUIRE(seeder.execute(ctx) == ActsExamples::ProcessCode::SUCCESS);
0132
0133   ActsExamples::ReadDataHandle<ActsExamples::MuonHoughMaxContainer>
0134       outputHandle{&seeder, "TestOutputHoughMaxima"};
0135   outputHandle.initialize(cfg.outHoughMax);
0136
0137   const ActsExamples::MuonHoughMaxContainer& maxima = outputHandle(ctx);
0138
0139   BOOST_REQUIRE_GT(maxima.size(), 0u);
0140
0141   bool foundExpectedMaximum = false;
0142
0143   double foundTanBeta = std::numeric_limits<double>::quiet_NaN();
0144   double foundInterceptY = std::numeric_limits<double>::quiet_NaN();
0145   for (const ActsExamples::MuonHoughMaximum& maximum : maxima) {
0146     const double tanTheta = maximum.tanBeta();
0147     const double interceptY = maximum.interceptY();
0148
0149     if (std::abs(tanTheta - trueTanTheta) < 0.02 &&
0150         std::abs(interceptY - trueInterceptY) < 20.0) {
0151       foundExpectedMaximum = true;
0152       foundTanBeta = tanTheta;
0153       foundInterceptY = interceptY;
0154       break;
0155     }
0156   }
0157
0158   /// Result acquired: -0.041608695652173948 -421.93644749999976
0159   BOOST_TEST_MESSAGE("Maximum coordinates (tanBeta, interceptY): "
0160                      << foundTanBeta << " " << foundInterceptY
0161                      << " truth: " << trueTanTheta << " " << trueInterceptY);
0162   BOOST_CHECK(foundExpectedMaximum);
0163 }
0164
0165 BOOST_AUTO_TEST_SUITE_END()
0166
0167 }  // namespace ActsTests