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