Back to home page

EIC code displayed by LXR
 
 

    

File indexing completed on 2025-01-18 09:12:47 

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

This page was automatically generated by the 2.3.7 LXR engine.

The LXR team
Valid CSS 2.1! Valid HTML 4.01!