Core/Propagator/MaterialCollectionTests.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/Direction.hpp"
0014 #include "Acts/Definitions/TrackParametrization.hpp"
0015 #include "Acts/Definitions/Units.hpp"
0016 #include "Acts/EventData/TrackParameters.hpp"
0017 #include "Acts/Geometry/GeometryContext.hpp"
0018 #include "Acts/Geometry/GeometryIdentifier.hpp"
0019 #include "Acts/MagneticField/ConstantBField.hpp"
0020 #include "Acts/MagneticField/MagneticFieldContext.hpp"
0021 #include "Acts/Propagator/ActorList.hpp"
0022 #include "Acts/Propagator/EigenStepper.hpp"
0023 #include "Acts/Propagator/MaterialInteractor.hpp"
0024 #include "Acts/Propagator/Navigator.hpp"
0025 #include "Acts/Propagator/Propagator.hpp"
0026 #include "Acts/Propagator/StraightLineStepper.hpp"
0027 #include "Acts/Surfaces/Surface.hpp"
0028 #include "ActsTests/CommonHelpers/CylindricalTrackingGeometry.hpp"
0029 #include "ActsTests/CommonHelpers/FloatComparisons.hpp"
0030
0031 #include <cmath>
0032 #include <iostream>
0033 #include <memory>
0034 #include <numbers>
0035 #include <optional>
0036 #include <random>
0037 #include <utility>
0038 #include <vector>
0039
0040 namespace bdata = boost::unit_test::data;
0041
0042 using namespace Acts;
0043 using namespace Acts::UnitLiterals;
0044
0045 namespace ActsTests {
0046
0047 // Create a test context
0048 GeometryContext tgContext = GeometryContext();
0049 MagneticFieldContext mfContext = MagneticFieldContext();
0050
0051 // Global definitions
0052 CylindricalTrackingGeometry cGeometry(tgContext);
0053 auto tGeometry = cGeometry();
0054
0055 using BField = ConstantBField;
0056 using EigenStepper = EigenStepper<>;
0057 using EigenPropagator = Propagator<EigenStepper, Navigator>;
0058 using StraightLinePropagator = Propagator<StraightLineStepper, Navigator>;
0059
0060 const double Bz = 2_T;
0061 auto bField = std::make_shared<BField>(Vector3{0, 0, Bz});
0062
0063 EigenStepper estepper(bField);
0064 Navigator esnavigator({tGeometry});
0065 EigenPropagator epropagator(std::move(estepper), std::move(esnavigator));
0066
0067 StraightLineStepper slstepper;
0068 Navigator slnavigator({tGeometry});
0069 StraightLinePropagator slpropagator(slstepper, std::move(slnavigator));
0070
0071 int ntests = 500;
0072 int skip = 0;
0073 bool debugMode = false;
0074
0075 /// the actual test nethod that runs the test can be used with several
0076 /// propagator types
0077 ///
0078 /// @tparam propagator_t is the actual propagator type
0079 ///
0080 /// @param prop is the propagator instance
0081 /// @param start the start parameters
0082 template <typename propagator_t>
0083 void runTest(const propagator_t& prop, const BoundTrackParameters& start) {
0084   // Action list and abort list
0085   using ActorList = ActorList<MaterialInteractor>;
0086
0087   using Options = typename propagator_t::template Options<ActorList>;
0088   Options fwdOptions(tgContext, mfContext);
0089   fwdOptions.stepping.maxStepSize = 25_cm;
0090   fwdOptions.pathLimit = 25_cm;
0091
0092   // get the material collector and configure it
0093   auto& fwdMaterialInteractor =
0094       fwdOptions.actorList.template get<MaterialInteractor>();
0095   fwdMaterialInteractor.recordInteractions = true;
0096   fwdMaterialInteractor.energyLoss = false;
0097   fwdMaterialInteractor.multipleScattering = false;
0098
0099   if (debugMode) {
0100     std::cout << ">>> Forward Propagation : start." << std::endl;
0101   }
0102   // forward material test
0103   const auto& fwdResult = prop.propagate(start, fwdOptions).value();
0104   const auto& fwdMaterial =
0105       fwdResult.template get<MaterialInteractor::result_type>();
0106   // check that the collected material is not zero
0107   BOOST_CHECK_NE(fwdMaterial.materialInX0, 0.);
0108   BOOST_CHECK_NE(fwdMaterial.materialInL0, 0.);
0109
0110   double fwdStepMaterialInX0 = 0.;
0111   double fwdStepMaterialInL0 = 0.;
0112   // check that the sum of all steps is the total material
0113   for (auto& mInteraction : fwdMaterial.materialInteractions) {
0114     fwdStepMaterialInX0 += mInteraction.materialSlab.thicknessInX0();
0115     fwdStepMaterialInL0 += mInteraction.materialSlab.thicknessInL0();
0116   }
0117   CHECK_CLOSE_REL(fwdMaterial.materialInX0, fwdStepMaterialInX0, 1e-3);
0118   CHECK_CLOSE_REL(fwdMaterial.materialInL0, fwdStepMaterialInL0, 1e-3);
0119
0120   // get the forward output to the screen
0121   if (debugMode) {
0122     // check if the surfaces are free
0123     std::cout << ">>> Material steps found on ..." << std::endl;
0124     for (auto& fwdStepsC : fwdMaterial.materialInteractions) {
0125       std::cout << "--> Surface with " << fwdStepsC.surface->geometryId()
0126                 << std::endl;
0127     }
0128   }
0129
0130   // backward material test
0131   Options bwdOptions(tgContext, mfContext);
0132   bwdOptions.stepping.maxStepSize = 25_cm;
0133   bwdOptions.pathLimit = -25_cm;
0134   bwdOptions.direction = Direction::Backward();
0135
0136   // get the material collector and configure it
0137   auto& bwdMaterialInteractor =
0138       bwdOptions.actorList.template get<MaterialInteractor>();
0139   bwdMaterialInteractor.recordInteractions = true;
0140   bwdMaterialInteractor.energyLoss = false;
0141   bwdMaterialInteractor.multipleScattering = false;
0142
0143   const auto& startSurface = start.referenceSurface();
0144
0145   if (debugMode) {
0146     std::cout << ">>> Backward Propagation : start." << std::endl;
0147   }
0148   const auto& bwdResult =
0149       prop.propagate(*fwdResult.endParameters, startSurface, bwdOptions)
0150           .value();
0151   if (debugMode) {
0152     std::cout << ">>> Backward Propagation : end." << std::endl;
0153   }
0154   const auto& bwdMaterial =
0155       bwdResult.template get<typename MaterialInteractor::result_type>();
0156   // check that the collected material is not zero
0157   BOOST_CHECK_NE(bwdMaterial.materialInX0, 0.);
0158   BOOST_CHECK_NE(bwdMaterial.materialInL0, 0.);
0159
0160   double bwdStepMaterialInX0 = 0.;
0161   double bwdStepMaterialInL0 = 0.;
0162   // check that the sum of all steps is the total material
0163   for (auto& mInteraction : bwdMaterial.materialInteractions) {
0164     bwdStepMaterialInX0 += mInteraction.materialSlab.thicknessInX0();
0165     bwdStepMaterialInL0 += mInteraction.materialSlab.thicknessInL0();
0166   }
0167   CHECK_CLOSE_REL(bwdMaterial.materialInX0, bwdStepMaterialInX0, 1e-3);
0168   CHECK_CLOSE_REL(bwdMaterial.materialInL0, bwdStepMaterialInL0, 1e-3);
0169
0170   // get the backward output to the screen
0171   if (debugMode) {
0172     // check if the surfaces are free
0173     std::cout << ">>> Material steps found on ..." << std::endl;
0174     for (auto& bwdStepsC : bwdMaterial.materialInteractions) {
0175       std::cout << "--> Surface with " << bwdStepsC.surface->geometryId()
0176                 << std::endl;
0177     }
0178   }
0179
0180   // forward-backward compatibility test
0181   BOOST_CHECK_EQUAL(bwdMaterial.materialInteractions.size(),
0182                     fwdMaterial.materialInteractions.size());
0183
0184   CHECK_CLOSE_REL(bwdMaterial.materialInX0, fwdMaterial.materialInX0, 1e-3);
0185   CHECK_CLOSE_REL(bwdMaterial.materialInL0, fwdMaterial.materialInL0, 1e-3);
0186
0187   // stepping from one surface to the next
0188   // now go from surface to surface and check
0189   Options fwdStepOptions(tgContext, mfContext);
0190   fwdStepOptions.stepping.maxStepSize = 25_cm;
0191   fwdStepOptions.pathLimit = 25_cm;
0192
0193   // get the material collector and configure it
0194   auto& fwdStepMaterialInteractor =
0195       fwdStepOptions.actorList.template get<MaterialInteractor>();
0196   fwdStepMaterialInteractor.recordInteractions = true;
0197   fwdStepMaterialInteractor.energyLoss = false;
0198   fwdStepMaterialInteractor.multipleScattering = false;
0199
0200   double fwdStepStepMaterialInX0 = 0.;
0201   double fwdStepStepMaterialInL0 = 0.;
0202
0203   if (debugMode) {
0204     // check if the surfaces are free
0205     std::cout << ">>> Forward steps to be processed sequentially ..."
0206               << std::endl;
0207     for (auto& fwdStepsC : fwdMaterial.materialInteractions) {
0208       std::cout << "--> Surface with " << fwdStepsC.surface->geometryId()
0209                 << std::endl;
0210     }
0211   }
0212
0213   // move forward step by step through the surfaces
0214   BoundTrackParameters sParameters = start;
0215   std::vector<BoundTrackParameters> stepParameters;
0216   for (auto& fwdSteps : fwdMaterial.materialInteractions) {
0217     if (debugMode) {
0218       std::cout << ">>> Forward step : "
0219                 << sParameters.referenceSurface().geometryId() << " --> "
0220                 << fwdSteps.surface->geometryId() << std::endl;
0221     }
0222
0223     // make a forward step
0224     const auto& fwdStep =
0225         prop.propagate(sParameters, (*fwdSteps.surface), fwdStepOptions)
0226             .value();
0227
0228     auto& fwdStepMaterial =
0229         fwdStep.template get<typename MaterialInteractor::result_type>();
0230     fwdStepStepMaterialInX0 += fwdStepMaterial.materialInX0;
0231     fwdStepStepMaterialInL0 += fwdStepMaterial.materialInL0;
0232
0233     if (fwdStep.endParameters.has_value()) {
0234       // make sure the parameters do not run out of scope
0235       stepParameters.push_back(*fwdStep.endParameters);
0236       sParameters = stepParameters.back();
0237     }
0238   }
0239   // final destination surface
0240   const Surface& dSurface = fwdResult.endParameters->referenceSurface();
0241
0242   if (debugMode) {
0243     std::cout << ">>> Forward step : "
0244               << sParameters.referenceSurface().geometryId() << " --> "
0245               << dSurface.geometryId() << std::endl;
0246   }
0247
0248   const auto& fwdStepFinal =
0249       prop.propagate(sParameters, dSurface, fwdStepOptions).value();
0250
0251   auto& fwdStepMaterial =
0252       fwdStepFinal.template get<typename MaterialInteractor::result_type>();
0253   fwdStepStepMaterialInX0 += fwdStepMaterial.materialInX0;
0254   fwdStepStepMaterialInL0 += fwdStepMaterial.materialInL0;
0255
0256   // forward-forward step compatibility test
0257   CHECK_CLOSE_REL(fwdStepStepMaterialInX0, fwdStepMaterialInX0, 1e-3);
0258   CHECK_CLOSE_REL(fwdStepStepMaterialInL0, fwdStepMaterialInL0, 1e-3);
0259
0260   // stepping from one surface to the next : backwards
0261   // now go from surface to surface and check
0262   Options bwdStepOptions(tgContext, mfContext);
0263   bwdStepOptions.stepping.maxStepSize = 25_cm;
0264   bwdStepOptions.pathLimit = -25_cm;
0265   bwdStepOptions.direction = Direction::Backward();
0266
0267   // get the material collector and configure it
0268   auto& bwdStepMaterialInteractor =
0269       bwdStepOptions.actorList.template get<MaterialInteractor>();
0270   bwdStepMaterialInteractor.recordInteractions = true;
0271   bwdStepMaterialInteractor.multipleScattering = false;
0272   bwdStepMaterialInteractor.energyLoss = false;
0273
0274   double bwdStepStepMaterialInX0 = 0.;
0275   double bwdStepStepMaterialInL0 = 0.;
0276
0277   if (debugMode) {
0278     // check if the surfaces are free
0279     std::cout << ">>> Backward steps to be processed sequentially ..."
0280               << std::endl;
0281     for (auto& bwdStepsC : bwdMaterial.materialInteractions) {
0282       std::cout << "--> Surface with " << bwdStepsC.surface->geometryId()
0283                 << std::endl;
0284     }
0285   }
0286
0287   // move forward step by step through the surfaces
0288   sParameters = *fwdResult.endParameters;
0289   for (auto& bwdSteps : bwdMaterial.materialInteractions) {
0290     if (debugMode) {
0291       std::cout << ">>> Backward step : "
0292                 << sParameters.referenceSurface().geometryId() << " --> "
0293                 << bwdSteps.surface->geometryId() << std::endl;
0294     }
0295     // make a forward step
0296     const auto& bwdStep =
0297         prop.propagate(sParameters, (*bwdSteps.surface), bwdStepOptions)
0298             .value();
0299
0300     auto& bwdStepMaterial =
0301         bwdStep.template get<typename MaterialInteractor::result_type>();
0302     bwdStepStepMaterialInX0 += bwdStepMaterial.materialInX0;
0303     bwdStepStepMaterialInL0 += bwdStepMaterial.materialInL0;
0304
0305     if (bwdStep.endParameters.has_value()) {
0306       // make sure the parameters do not run out of scope
0307       stepParameters.push_back(*bwdStep.endParameters);
0308       sParameters = stepParameters.back();
0309     }
0310   }
0311   // final destination surface
0312   const Surface& dbSurface = start.referenceSurface();
0313
0314   if (debugMode) {
0315     std::cout << ">>> Backward step : "
0316               << sParameters.referenceSurface().geometryId() << " --> "
0317               << dSurface.geometryId() << std::endl;
0318   }
0319
0320   const auto& bwdStepFinal =
0321       prop.propagate(sParameters, dbSurface, bwdStepOptions).value();
0322
0323   auto& bwdStepMaterial =
0324       bwdStepFinal.template get<typename MaterialInteractor::result_type>();
0325   bwdStepStepMaterialInX0 += bwdStepMaterial.materialInX0;
0326   bwdStepStepMaterialInL0 += bwdStepMaterial.materialInL0;
0327
0328   // backward-backward step compatibility test
0329   CHECK_CLOSE_REL(bwdStepStepMaterialInX0, bwdStepMaterialInX0, 1e-3);
0330   CHECK_CLOSE_REL(bwdStepStepMaterialInL0, bwdStepMaterialInL0, 1e-3);
0331
0332   // Test the material affects the covariance into the right direction
0333   // get the material collector and configure it
0334   auto& covfwdMaterialInteractor =
0335       fwdOptions.actorList.template get<MaterialInteractor>();
0336   covfwdMaterialInteractor.recordInteractions = false;
0337   covfwdMaterialInteractor.energyLoss = true;
0338   covfwdMaterialInteractor.multipleScattering = true;
0339
0340   // forward material test
0341   const auto& covfwdResult = prop.propagate(start, fwdOptions).value();
0342
0343   BOOST_CHECK_LE(
0344       start.covariance()->determinant(),
0345       covfwdResult.endParameters->covariance().value().determinant());
0346 }
0347
0348 BOOST_AUTO_TEST_SUITE(PropagatorSuite)
0349
0350 // This test case checks that no segmentation fault appears
0351 // - this tests the collection of surfaces
0352 BOOST_DATA_TEST_CASE(
0353     test_material_collector,
0354     bdata::random((bdata::engine = std::mt19937(), bdata::seed = 20,
0355                    bdata::distribution = std::uniform_real_distribution<double>(
0356                        0.5_GeV, 10_GeV))) ^
0357         bdata::random(
0358             (bdata::engine = std::mt19937(), bdata::seed = 21,
0359              bdata::distribution = std::uniform_real_distribution<double>(
0360                  -std::numbers::pi, std::numbers::pi))) ^
0361         bdata::random(
0362             (bdata::engine = std::mt19937(), bdata::seed = 22,
0363              bdata::distribution = std::uniform_real_distribution<double>(
0364                  1., std::numbers::pi - 1.))) ^
0365         bdata::random((bdata::engine = std::mt19937(), bdata::seed = 23,
0366                        bdata::distribution =
0367                            std::uniform_int_distribution<std::uint8_t>(0, 1))) ^
0368         bdata::xrange(ntests),
0369     pT, phi, theta, charge, index) {
0370   if (index < skip) {
0371     return;
0372   }
0373
0374   double p = pT / std::sin(theta);
0375   double q = -1 + 2 * charge;
0376
0377   // define start parameters
0378   BoundSquareMatrix cov;
0379   // take some major correlations (off-diagonals)
0380   // clang-format off
0381     cov <<
0382      10_mm, 0, 0.123, 0, 0.5, 0,
0383      0, 10_mm, 0, 0.162, 0, 0,
0384      0.123, 0, 0.1, 0, 0, 0,
0385      0, 0.162, 0, 0.1, 0, 0,
0386      0.5, 0, 0, 0, 1_e / 10_GeV, 0,
0387      0, 0, 0, 0, 0, 1_us;
0388   // clang-format on
0389   BoundTrackParameters start = BoundTrackParameters::createCurvilinear(
0390       Vector4::Zero(), phi, theta, q / p, cov, ParticleHypothesis::pion());
0391
0392   runTest(epropagator, start);
0393   runTest(slpropagator, start);
0394 }
0395
0396 BOOST_AUTO_TEST_SUITE_END()
0397
0398 }  // namespace ActsTests