Core/SpacePointFormation/SpacePointBuilderTests.cpp

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