Core/Propagator/MaterialCollectionTests.cpp

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