Core/SpacePointFormation/SpacePointBuilderTests.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/data/test_case.hpp>
0010 #include <boost/test/unit_test.hpp>
0011
0012 #include "Acts/Definitions/Algebra.hpp"
0013 #include "Acts/Definitions/TrackParametrization.hpp"
0014 #include "Acts/Definitions/Units.hpp"
0015 #include "Acts/EventData/SourceLink.hpp"
0016 #include "Acts/EventData/TrackParameters.hpp"
0017 #include "Acts/EventData/detail/TestSourceLink.hpp"
0018 #include "Acts/Geometry/GeometryContext.hpp"
0019 #include "Acts/Geometry/GeometryIdentifier.hpp"
0020 #include "Acts/Geometry/TrackingGeometry.hpp"
0021 #include "Acts/MagneticField/ConstantBField.hpp"
0022 #include "Acts/MagneticField/MagneticFieldContext.hpp"
0023 #include "Acts/Propagator/EigenStepper.hpp"
0024 #include "Acts/Propagator/Navigator.hpp"
0025 #include "Acts/Propagator/Propagator.hpp"
0026 #include "Acts/SpacePointFormation/SpacePointBuilder.hpp"
0027 #include "Acts/SpacePointFormation/SpacePointBuilderConfig.hpp"
0028 #include "Acts/SpacePointFormation/SpacePointBuilderOptions.hpp"
0029 #include "Acts/Surfaces/Surface.hpp"
0030 #include "Acts/Tests/CommonHelpers/CubicTrackingGeometry.hpp"
0031 #include "Acts/Tests/CommonHelpers/MeasurementsCreator.hpp"
0032 #include "Acts/Tests/CommonHelpers/TestSpacePoint.hpp"
0033
0034 #include <iostream>
0035 #include <iterator>
0036 #include <memory>
0037 #include <optional>
0038 #include <random>
0039 #include <utility>
0040 #include <vector>
0041
0042 namespace bdata = boost::unit_test::data;
0043 using namespace Acts::UnitLiterals;
0044
0045 namespace Acts::Test {
0046
0047 using TestSourceLink = detail::Test::TestSourceLink;
0048 using ConstantFieldStepper = EigenStepper<>;
0049 using ConstantFieldPropagator = Propagator<ConstantFieldStepper, Navigator>;
0050
0051 /// Construct initial track parameters.
0052 CurvilinearTrackParameters makeParameters(double phi, double theta, double p,
0053                                           double q) {
0054   // create covariance matrix from reasonable standard deviations
0055   BoundVector stddev;
0056   stddev[eBoundLoc0] = 100_um;
0057   stddev[eBoundLoc1] = 100_um;
0058   stddev[eBoundTime] = 25_ns;
0059   stddev[eBoundPhi] = 2_degree;
0060   stddev[eBoundTheta] = 2_degree;
0061   stddev[eBoundQOverP] = 1 / 100_GeV;
0062   BoundSquareMatrix cov = stddev.cwiseProduct(stddev).asDiagonal();
0063   // Let the particle start from the origin
0064   Vector4 mPos4(-3_m, 0., 0., 0.);
0065   return CurvilinearTrackParameters(mPos4, phi, theta, q / p, cov,
0066                                     ParticleHypothesis::pionLike(q));
0067 }
0068
0069 std::pair<Vector3, Vector3> stripEnds(
0070     const std::shared_ptr<const TrackingGeometry>& geo,
0071     const GeometryContext& gctx, const SourceLink& slink,
0072     const double stripFrac = 0.4) {
0073   auto testslink = slink.get<TestSourceLink>();
0074   const auto lpos = testslink.parameters;
0075
0076   Vector3 globalFakeMom(1, 1, 1);
0077   const auto geoId = testslink.m_geometryId;
0078   const Surface* surface = geo->findSurface(geoId);
0079
0080   const double stripLength = 40.;
0081   const double end1x = lpos[0] + stripLength * stripFrac;
0082   const double end1y = lpos[1];
0083   const double end2x = lpos[0] - stripLength * (1 - stripFrac);
0084   const double end2y = lpos[1];
0085   const Vector2 lpos1(end1x, end1y);
0086   const Vector2 lpos2(end2x, end2y);
0087
0088   auto gPos1 = surface->localToGlobal(gctx, lpos1, globalFakeMom);
0089   auto gPos2 = surface->localToGlobal(gctx, lpos2, globalFakeMom);
0090
0091   return {gPos1, gPos2};
0092 }
0093
0094 // Create a test context
0095 GeometryContext tgContext = GeometryContext();
0096
0097 const GeometryContext geoCtx;
0098 const MagneticFieldContext magCtx;
0099
0100 // detector geometry
0101 CubicTrackingGeometry geometryStore(geoCtx);
0102 const auto geometry = geometryStore();
0103
0104 // detector resolutions
0105 const MeasurementResolution resPixel = {MeasurementType::eLoc01,
0106                                         {25_um, 50_um}};
0107 const MeasurementResolution resStrip = {MeasurementType::eLoc01,
0108                                         {100_um, 100_um}};
0109 const MeasurementResolutionMap resolutions = {
0110     {GeometryIdentifier().setVolume(2), resPixel},
0111     {GeometryIdentifier().setVolume(3).setLayer(2), resStrip},
0112     {GeometryIdentifier().setVolume(3).setLayer(4), resStrip},
0113     {GeometryIdentifier().setVolume(3).setLayer(6), resStrip},
0114     {GeometryIdentifier().setVolume(3).setLayer(8), resStrip},
0115 };
0116
0117 std::default_random_engine rng(42);
0118
0119 BOOST_DATA_TEST_CASE(SpacePointBuilder_basic, bdata::xrange(1), index) {
0120   (void)index;
0121
0122   double phi = 5._degree;
0123   double theta = 95._degree;
0124   double p = 50._GeV;
0125   double q = 1;
0126
0127   Navigator navigator({
0128       geometry,
0129       true,  // sensitive
0130       true,  // material
0131       false  // passive
0132   });
0133   auto field = std::make_shared<ConstantBField>(Vector3(0.0, 0.0, 2._T));
0134   ConstantFieldStepper stepper(std::move(field));
0135
0136   ConstantFieldPropagator propagator(std::move(stepper), std::move(navigator));
0137   auto start = makeParameters(phi, theta, p, q);
0138
0139   auto measurements =
0140       createMeasurements(propagator, geoCtx, magCtx, start, resolutions, rng);
0141
0142   const auto sourceLinks = measurements.sourceLinks;
0143
0144   std::vector<SourceLink> frontSourceLinks;
0145   std::vector<SourceLink> backSourceLinks;
0146   std::vector<SourceLink> singleHitSourceLinks;
0147
0148   std::vector<const Vector3*> frontStripEnds;
0149   std::vector<const Vector3*> backStripEnds;
0150
0151   for (auto& sl : sourceLinks) {
0152     const auto geoId = sl.m_geometryId;
0153     const auto volumeId = geoId.volume();
0154     if (volumeId == 2) {  // pixel type detector
0155       singleHitSourceLinks.emplace_back(SourceLink{sl});
0156     } else if (volumeId == 3) {  // strip type detector
0157
0158       const auto layerId = geoId.layer();
0159
0160       if (layerId == 2 || layerId == 6) {
0161         frontSourceLinks.emplace_back(SourceLink{sl});
0162       } else if (layerId == 4 || layerId == 8) {
0163         backSourceLinks.emplace_back(SourceLink{sl});
0164       }
0165     }  // volume 3 (strip detector)
0166   }
0167
0168   BOOST_CHECK_EQUAL(frontSourceLinks.size(), 2);
0169   BOOST_CHECK_EQUAL(backSourceLinks.size(), 2);
0170
0171   Vector3 vertex = Vector3(-3_m, 0., 0.);
0172
0173   auto spConstructor = [](const Vector3& pos, const std::optional<double>& t,
0174                           const Vector2& cov, const std::optional<double>& covT,
0175                           boost::container::static_vector<SourceLink, 2> slinks)
0176       -> TestSpacePoint {
0177     return TestSpacePoint(pos, t, cov[0], cov[1], covT, std::move(slinks));
0178   };
0179
0180   auto spBuilderConfig = SpacePointBuilderConfig();
0181   spBuilderConfig.trackingGeometry = geometry;
0182
0183   TestSourceLink::SurfaceAccessor surfaceAccessor{*geometry};
0184   spBuilderConfig.slSurfaceAccessor
0185       .connect<&TestSourceLink::SurfaceAccessor::operator()>(&surfaceAccessor);
0186
0187   auto spBuilder =
0188       SpacePointBuilder<TestSpacePoint>(spBuilderConfig, spConstructor);
0189
0190   // for cosmic  without vertex constraint, usePerpProj = true
0191   auto spBuilderConfig_perp = SpacePointBuilderConfig();
0192   spBuilderConfig_perp.trackingGeometry = geometry;
0193   spBuilderConfig_perp.slSurfaceAccessor
0194       .connect<&TestSourceLink::SurfaceAccessor::operator()>(&surfaceAccessor);
0195
0196   spBuilderConfig_perp.usePerpProj = true;
0197
0198   auto spBuilder_perp =
0199       SpacePointBuilder<TestSpacePoint>(spBuilderConfig_perp, spConstructor);
0200
0201   TestSpacePointContainer spacePoints;
0202   TestSpacePointContainer spacePoints_extra;
0203
0204   auto accessor = [](const SourceLink& slink) {
0205     auto testslink = slink.get<TestSourceLink>();
0206     BoundVector param;
0207     param.setZero();
0208     param[eBoundLoc0] = testslink.parameters[eBoundLoc0];
0209     param[eBoundLoc1] = testslink.parameters[eBoundLoc1];
0210
0211     BoundSquareMatrix cov = BoundSquareMatrix::Zero();
0212     cov.topLeftCorner<2, 2>() = testslink.covariance;
0213
0214     return std::make_pair(param, cov);
0215   };
0216
0217   for (auto& sl : singleHitSourceLinks) {
0218     std::vector<SourceLink> slinks;
0219     slinks.emplace_back(sl);
0220     SpacePointBuilderOptions spOpt;
0221     spOpt.vertex = vertex;
0222     spOpt.paramCovAccessor = accessor;
0223     spBuilder.buildSpacePoint(geoCtx, slinks, spOpt,
0224                               std::back_inserter(spacePoints));
0225   }
0226   BOOST_CHECK_EQUAL(spacePoints.size(), 2);
0227   std::vector<std::pair<SourceLink, SourceLink>> slinkPairs;
0228
0229   // strip SP building
0230
0231   StripPairOptions pairOpt;
0232   pairOpt.paramCovAccessor = accessor;
0233
0234   spBuilder.makeSourceLinkPairs(tgContext, frontSourceLinks, backSourceLinks,
0235                                 slinkPairs, pairOpt);
0236
0237   BOOST_CHECK_EQUAL(slinkPairs.size(), 2);
0238
0239   for (auto& slinkPair : slinkPairs) {
0240     const std::pair<Vector3, Vector3> end1 =
0241         stripEnds(geometry, geoCtx, slinkPair.first);
0242     const std::pair<Vector3, Vector3> end2 =
0243         stripEnds(geometry, geoCtx, slinkPair.second);
0244
0245     std::shared_ptr<const TestSpacePoint> spacePoint = nullptr;
0246
0247     auto strippair = std::make_pair(end1, end2);
0248     std::vector<SourceLink> slinks;
0249     slinks.emplace_back(slinkPair.first);
0250     slinks.emplace_back(slinkPair.second);
0251
0252     SpacePointBuilderOptions spOpt{strippair, accessor};
0253
0254     // nominal strip sp building
0255     spBuilder.buildSpacePoint(geoCtx, slinks, spOpt,
0256                               std::back_inserter(spacePoints));
0257
0258     // sp building without vertex constraint
0259     spBuilder_perp.buildSpacePoint(geoCtx, slinks, spOpt,
0260                                    std::back_inserter(spacePoints));
0261
0262     // put measurements slightly outside strips to test recovery
0263     const std::pair<Vector3, Vector3> end3 =
0264         stripEnds(geometry, geoCtx, slinkPair.first, 1.01);
0265     const std::pair<Vector3, Vector3> end4 =
0266         stripEnds(geometry, geoCtx, slinkPair.second, 1.02);
0267     // the other side of the strips
0268     const std::pair<Vector3, Vector3> end5 =
0269         stripEnds(geometry, geoCtx, slinkPair.first, -0.01);
0270     const std::pair<Vector3, Vector3> end6 =
0271         stripEnds(geometry, geoCtx, slinkPair.second, -0.02);
0272
0273     auto spBuilderConfig_badStrips = SpacePointBuilderConfig();
0274
0275     spBuilderConfig_badStrips.trackingGeometry = geometry;
0276     spBuilderConfig_badStrips.slSurfaceAccessor
0277         .connect<&TestSourceLink::SurfaceAccessor::operator()>(
0278             &surfaceAccessor);
0279
0280     auto spBuilder_badStrips = SpacePointBuilder<TestSpacePoint>(
0281         spBuilderConfig_badStrips, spConstructor);
0282     // sp building with the recovery method
0283     SpacePointBuilderOptions spOpt_badStrips1{std::make_pair(end3, end4),
0284                                               accessor};
0285     spOpt_badStrips1.vertex = vertex;
0286     spOpt_badStrips1.stripLengthTolerance = 0.0001;
0287     spOpt_badStrips1.stripLengthGapTolerance = 50.;
0288     spBuilder_badStrips.buildSpacePoint(geoCtx, slinks, spOpt_badStrips1,
0289                                         std::back_inserter(spacePoints_extra));
0290
0291     SpacePointBuilderOptions spOpt_badStrips2{std::make_pair(end5, end6),
0292                                               accessor};
0293     spOpt_badStrips2.vertex = vertex;
0294     spOpt_badStrips2.stripLengthTolerance = 0.0001;
0295     spOpt_badStrips2.stripLengthGapTolerance = 50.;
0296     spBuilder_badStrips.buildSpacePoint(geoCtx, slinks, spOpt_badStrips2,
0297                                         std::back_inserter(spacePoints_extra));
0298   }
0299
0300   for (auto& sp : spacePoints) {
0301     std::cout << "space point (" << sp.x() << " " << sp.y() << " " << sp.z()
0302               << ") var (r,z): " << sp.varianceR() << " " << sp.varianceZ()
0303               << std::endl;
0304   }
0305   std::cout << "space points produced with bad strips:" << std::endl;
0306   for (auto& sp : spacePoints_extra) {
0307     std::cout << "space point (" << sp.x() << " " << sp.y() << " " << sp.z()
0308               << ") var (r,z): " << sp.varianceR() << " " << sp.varianceZ()
0309               << std::endl;
0310   }
0311
0312   BOOST_CHECK_EQUAL(spacePoints.size(), 6);
0313 }
0314
0315 }  // namespace Acts::Test