Core/Visualization/EventDataView3DBase.hpp

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 #pragma once
0010
0011 #include "Acts/Definitions/Algebra.hpp"
0012 #include "Acts/EventData/SourceLink.hpp"
0013 #include "Acts/EventData/VectorMultiTrajectory.hpp"
0014 #include "Acts/EventData/VectorTrackContainer.hpp"
0015 #include "Acts/EventData/detail/TestSourceLink.hpp"
0016 #include "Acts/Geometry/CuboidVolumeBuilder.hpp"
0017 #include "Acts/Geometry/GeometryContext.hpp"
0018 #include "Acts/Geometry/TrackingGeometry.hpp"
0019 #include "Acts/Geometry/TrackingGeometryBuilder.hpp"
0020 #include "Acts/MagneticField/ConstantBField.hpp"
0021 #include "Acts/MagneticField/MagneticFieldContext.hpp"
0022 #include "Acts/Material/HomogeneousSurfaceMaterial.hpp"
0023 #include "Acts/Propagator/EigenStepper.hpp"
0024 #include "Acts/Propagator/Navigator.hpp"
0025 #include "Acts/Propagator/Propagator.hpp"
0026 #include "Acts/Surfaces/PlaneSurface.hpp"
0027 #include "Acts/Surfaces/RectangleBounds.hpp"
0028 #include "Acts/TrackFitting/GainMatrixSmoother.hpp"
0029 #include "Acts/TrackFitting/GainMatrixUpdater.hpp"
0030 #include "Acts/TrackFitting/KalmanFitter.hpp"
0031 #include "Acts/Utilities/CalibrationContext.hpp"
0032 #include "Acts/Visualization/EventDataView3D.hpp"
0033 #include "Acts/Visualization/IVisualization3D.hpp"
0034 #include "ActsTests/CommonHelpers/DetectorElementStub.hpp"
0035 #include "ActsTests/CommonHelpers/PredefinedMaterials.hpp"
0036
0037 #include <cmath>
0038 #include <optional>
0039 #include <random>
0040 #include <sstream>
0041 #include <string>
0042
0043 using Acts::VectorHelpers::makeVector4;
0044
0045 namespace Acts::EventDataView3DTest {
0046
0047 using Covariance = BoundMatrix;
0048
0049 std::normal_distribution<double> gauss(0., 1.);
0050 std::default_random_engine generator(42);
0051
0052 /// A function that creates a simple telescope detector with surfaces for the
0053 /// EventDataView3D tests
0054 ///
0055 /// @param helper The visualization helper
0056 /// @param surfaces Reference to the surfaces, because we need them outside
0057 /// @param detector The tracking geometry that will be filled
0058 /// @param nSurfaces Number of surfaces to generate
0059 ///
0060 /// @return an overall string including all written output
0061 void createDetector(GeometryContext& tgContext,
0062                     std::vector<const Surface*>& surfaces,
0063                     std::shared_ptr<const TrackingGeometry>& detector,
0064                     const std::size_t nSurfaces = 7) {
0065   using namespace UnitLiterals;
0066
0067   if (nSurfaces < 1) {
0068     throw std::invalid_argument("At least 1 surfaces needs to be created.");
0069   }
0070
0071   // Construct the rotation
0072   RotationMatrix3 rotation = RotationMatrix3::Identity();
0073   double rotationAngle = 90_degree;
0074   Vector3 xPos(std::cos(rotationAngle), 0., std::sin(rotationAngle));
0075   Vector3 yPos(0., 1., 0.);
0076   Vector3 zPos(-std::sin(rotationAngle), 0., std::cos(rotationAngle));
0077   rotation.col(0) = xPos;
0078   rotation.col(1) = yPos;
0079   rotation.col(2) = zPos;
0080
0081   // Boundaries of the surfaces
0082   const auto rBounds =
0083       std::make_shared<const RectangleBounds>(RectangleBounds(50_mm, 50_mm));
0084
0085   // Material of the surfaces
0086   MaterialSlab matProp(ActsTests::makeSilicon(), 0.5_mm);
0087   const auto surfaceMaterial =
0088       std::make_shared<HomogeneousSurfaceMaterial>(matProp);
0089
0090   // Set translation vectors
0091   std::vector<Vector3> translations;
0092   translations.reserve(nSurfaces);
0093   for (unsigned int i = 0; i < nSurfaces; i++) {
0094     translations.push_back({i * 100_mm - 300_mm, 0., 0.});
0095   }
0096
0097   // Construct layer configs
0098   std::vector<CuboidVolumeBuilder::LayerConfig> lConfs;
0099   lConfs.reserve(nSurfaces);
0100   for (unsigned int i = 0; i < translations.size(); i++) {
0101     CuboidVolumeBuilder::SurfaceConfig sConf;
0102     sConf.position = translations[i];
0103     sConf.rotation = rotation;
0104     sConf.rBounds = rBounds;
0105     sConf.surMat = surfaceMaterial;
0106     // The thickness to construct the associated detector element
0107     sConf.thickness = 1._um;
0108     sConf.detElementConstructor =
0109         [](const Transform3& trans,
0110            const std::shared_ptr<const RectangleBounds>& bounds,
0111            double thickness) {
0112           return new ActsTests::DetectorElementStub(trans, bounds, thickness);
0113         };
0114     CuboidVolumeBuilder::LayerConfig lConf;
0115     lConf.surfaceCfg = {sConf};
0116     lConfs.push_back(lConf);
0117   }
0118
0119   // Construct volume config
0120   CuboidVolumeBuilder::VolumeConfig vConf;
0121   vConf.position = {0., 0., 0.};
0122   vConf.length = {1.2_m, 1._m, 1._m};
0123   vConf.layerCfg = lConfs;
0124   vConf.name = "Tracker";
0125
0126   // Construct volume builder config
0127   CuboidVolumeBuilder::Config conf;
0128   conf.position = {0., 0., 0.};
0129   conf.length = {1.2_m, 1._m, 1._m};
0130   conf.volumeCfg = {vConf};  // one volume
0131
0132   // Build detector
0133   std::cout << "Build the detector" << std::endl;
0134   CuboidVolumeBuilder cvb;
0135   cvb.setConfig(conf);
0136   TrackingGeometryBuilder::Config tgbCfg;
0137   tgbCfg.trackingVolumeBuilders.push_back(
0138       [=](const auto& context, const auto& inner, const auto& vb) {
0139         return cvb.trackingVolume(context, inner, vb);
0140       });
0141   TrackingGeometryBuilder tgb(tgbCfg);
0142   detector = tgb.trackingGeometry(tgContext);
0143
0144   // Get the surfaces;
0145   surfaces.reserve(nSurfaces);
0146   detector->visitSurfaces([&](const Surface* surface) {
0147     if (surface != nullptr && surface->isSensitive()) {
0148       std::cout << "surface " << surface->geometryId() << " placed at: ("
0149                 << surface->center(tgContext).transpose() << " )" << std::endl;
0150       surfaces.push_back(surface);
0151     }
0152   });
0153   std::cout << "There are " << surfaces.size() << " surfaces" << std::endl;
0154 }
0155
0156 /// Helper method to visualize all types of surfaces
0157 ///
0158 /// @param helper The visualization helper
0159 ///
0160 /// @return an overall string including all written output
0161 static inline std::string testBoundTrackParameters(IVisualization3D& helper) {
0162   std::stringstream ss;
0163
0164   ViewConfig pcolor{.color = {20, 120, 20}};
0165   ViewConfig scolor{.color = {235, 198, 52}};
0166
0167   auto gctx = GeometryContext::dangerouslyDefaultConstruct();
0168   auto identity = Transform3::Identity();
0169
0170   // rectangle and plane
0171   auto rectangle = std::make_shared<RectangleBounds>(15., 15.);
0172   auto plane = Surface::makeShared<PlaneSurface>(identity, rectangle);
0173
0174   double momentumScale = 0.005;
0175   double localErrorScale = 10.;
0176   double directionErrorScale = 1000.;
0177
0178   // now create parameters on this surface
0179   // l_x, l_y, phi, theta, q/p (1/p), t
0180   std::array<double, 6> pars_array = {
0181       {-0.1234, 4.8765, 0.45, 0.128, 0.001, 21.}};
0182
0183   BoundVector pars = BoundVector::Zero();
0184   pars << pars_array[0], pars_array[1], pars_array[2], pars_array[3],
0185       pars_array[4], pars_array[5];
0186
0187   Covariance cov = Covariance::Zero();
0188   cov << 0.25, 0.0042, -0.00076, 6.156e-06, -2.11e-07, 0, 0.0042, 0.859,
0189       -0.000173, 0.000916, -4.017e-08, 0, -0.00076, -0.000173, 2.36e-04,
0190       -2.76e-07, 1.12e-08, 0, 6.15e-06, 0.000916, -2.76e-07, 8.84e-04,
0191       -2.85e-11, 0, -2.11 - 07, -4.017e-08, 1.123e-08, -2.85 - 11, 1.26e-10, 0,
0192       0, 0, 0, 0, 0, 1;
0193
0194   EventDataView3D::drawBoundTrackParameters(
0195       helper,
0196       BoundTrackParameters(plane, pars, std::move(cov),
0197                            ParticleHypothesis::pion()),
0198       gctx, momentumScale, localErrorScale, directionErrorScale, pcolor,
0199       scolor);
0200
0201   helper.write("EventData_BoundAtPlaneParameters");
0202   helper.write(ss);
0203
0204   return ss.str();
0205 }
0206
0207 /// Helper method to visualize measurements
0208 ///
0209 /// @param helper The visualization helper
0210 ///
0211 /// @return an overall string including all written output
0212 static inline std::string testMeasurement(IVisualization3D& helper,
0213                                           const double localErrorScale = 100.) {
0214   using namespace UnitLiterals;
0215   std::stringstream ss;
0216
0217   // Create a test context
0218   GeometryContext tgContext = GeometryContext::dangerouslyDefaultConstruct();
0219
0220   // Create a detector
0221   const std::size_t nSurfaces = 7;
0222   std::vector<const Surface*> surfaces;
0223   std::shared_ptr<const TrackingGeometry> detector;
0224   createDetector(tgContext, surfaces, detector, nSurfaces);
0225
0226   // Create measurements (assuming they are for a linear track parallel to
0227   // global x-axis)
0228   std::cout << "Creating measurements:" << std::endl;
0229   std::vector<detail::Test::TestSourceLink> sourceLinks;
0230   sourceLinks.reserve(nSurfaces);
0231   Vector2 lPosCenter{5_mm, 5_mm};
0232   Vector2 resolution{200_um, 150_um};
0233   SquareMatrix2 cov2D = resolution.cwiseProduct(resolution).asDiagonal();
0234   for (const auto& surface : surfaces) {
0235     // 2D measurements
0236     Vector2 loc = lPosCenter;
0237     loc[0] += resolution[0] * gauss(generator);
0238     loc[1] += resolution[1] * gauss(generator);
0239     sourceLinks.emplace_back(detail::Test::TestSourceLink{
0240         eBoundLoc0, eBoundLoc1, loc, cov2D, surface->geometryId()});
0241   }
0242
0243   ViewConfig mcolor{.color = {255, 145, 48}};
0244   mcolor.offset = 0.01;
0245
0246   // Draw the measurements
0247   std::cout << "Draw the measurements" << std::endl;
0248   //  auto singleMeasurement = sourceLinks[0];
0249   for (auto& singleMeasurement : sourceLinks) {
0250     auto cov = singleMeasurement.covariance;
0251     auto lposition = singleMeasurement.parameters;
0252
0253     auto surf = detector->findSurface(singleMeasurement.m_geometryId);
0254     auto transf = surf->localToGlobalTransform(tgContext);
0255
0256     EventDataView3D::drawMeasurement(helper, lposition, cov, transf,
0257                                      localErrorScale, mcolor);
0258   }
0259
0260   helper.write("EventData_Measurement");
0261   helper.write(ss);
0262
0263   return ss.str();
0264 }
0265
0266 /// Helper method to visualize a MultiTrajectory
0267 ///
0268 /// @param helper The visualization helper
0269 ///
0270 /// @return an overall string including all written output
0271 static inline std::string testMultiTrajectory(IVisualization3D& helper) {
0272   using namespace UnitLiterals;
0273   std::stringstream ss;
0274
0275   // Create a test context
0276   GeometryContext tgContext = GeometryContext::dangerouslyDefaultConstruct();
0277   MagneticFieldContext mfContext = MagneticFieldContext();
0278   CalibrationContext calContext = CalibrationContext();
0279
0280   // Create a detector
0281   const std::size_t nSurfaces = 7;
0282   std::vector<const Surface*> surfaces;
0283   std::shared_ptr<const TrackingGeometry> detector;
0284   createDetector(tgContext, surfaces, detector, nSurfaces);
0285
0286   // Create measurements (assuming they are for a linear track parallel to
0287   // global x-axis)
0288   std::cout << "Creating measurements:" << std::endl;
0289   std::vector<Acts::SourceLink> sourceLinks;
0290   sourceLinks.reserve(nSurfaces);
0291   Vector2 lPosCenter{5_mm, 5_mm};
0292   Vector2 resolution{200_um, 150_um};
0293   SquareMatrix2 cov2D = resolution.cwiseProduct(resolution).asDiagonal();
0294   for (const auto& surface : surfaces) {
0295     // 2D measurements
0296     Vector2 loc = lPosCenter;
0297     loc[0] += resolution[0] * gauss(generator);
0298     loc[1] += resolution[1] * gauss(generator);
0299     sourceLinks.emplace_back(detail::Test::TestSourceLink{
0300         eBoundLoc0, eBoundLoc1, loc, cov2D, surface->geometryId()});
0301   }
0302
0303   // The KalmanFitter - we use the eigen stepper for covariance transport
0304   std::cout << "Construct KalmanFitter and perform fit" << std::endl;
0305   Navigator::Config cfg{detector};
0306   cfg.resolvePassive = false;
0307   cfg.resolveMaterial = true;
0308   cfg.resolveSensitive = true;
0309   Navigator rNavigator(cfg);
0310
0311   // Configure propagation with deactivated B-field
0312   auto bField = std::make_shared<ConstantBField>(Vector3(0., 0., 0.));
0313   using RecoStepper = EigenStepper<>;
0314   RecoStepper rStepper(bField);
0315   using RecoPropagator = Propagator<RecoStepper, Navigator>;
0316   RecoPropagator rPropagator(rStepper, rNavigator);
0317
0318   // Set initial parameters for the particle track
0319   Covariance cov;
0320   cov << std::pow(100_um, 2), 0., 0., 0., 0., 0., 0., std::pow(100_um, 2), 0.,
0321       0., 0., 0., 0., 0., 0.0025, 0., 0., 0., 0., 0., 0., 0.0025, 0., 0., 0.,
0322       0., 0., 0., 0.01, 0., 0., 0., 0., 0., 0., 1.;
0323   Vector3 rPos(-350._mm, 100_um * gauss(generator), 100_um * gauss(generator));
0324   Vector3 rDir(1, 0.025 * gauss(generator), 0.025 * gauss(generator));
0325   BoundTrackParameters rStart = BoundTrackParameters::createCurvilinear(
0326       makeVector4(rPos, 42_ns), rDir, 1_e / 1_GeV, cov,
0327       ParticleHypothesis::pion());
0328
0329   const Surface* rSurface = &rStart.referenceSurface();
0330
0331   using KalmanFitter = KalmanFitter<RecoPropagator, VectorMultiTrajectory>;
0332
0333   KalmanFitter kFitter(rPropagator);
0334
0335   auto logger = getDefaultLogger("KalmanFilter", Logging::WARNING);
0336
0337   Acts::GainMatrixUpdater kfUpdater;
0338   Acts::GainMatrixSmoother kfSmoother;
0339
0340   KalmanFitterExtensions<VectorMultiTrajectory> extensions;
0341   extensions.calibrator.connect<
0342       &detail::Test::testSourceLinkCalibrator<VectorMultiTrajectory>>();
0343   extensions.updater
0344       .connect<&Acts::GainMatrixUpdater::operator()<VectorMultiTrajectory>>(
0345           &kfUpdater);
0346   extensions.smoother
0347       .connect<&Acts::GainMatrixSmoother::operator()<VectorMultiTrajectory>>(
0348           &kfSmoother);
0349
0350   detail::Test::TestSourceLink::SurfaceAccessor surfaceAccessor{*detector};
0351   extensions.surfaceAccessor
0352       .connect<&detail::Test::TestSourceLink::SurfaceAccessor::operator()>(
0353           &surfaceAccessor);
0354
0355   KalmanFitterOptions kfOptions(tgContext, mfContext, calContext, extensions,
0356                                 PropagatorPlainOptions(tgContext, mfContext),
0357                                 rSurface);
0358
0359   Acts::TrackContainer tracks{Acts::VectorTrackContainer{},
0360                               Acts::VectorMultiTrajectory{}};
0361
0362   // Fit the track
0363   auto fitRes = kFitter.fit(sourceLinks.begin(), sourceLinks.end(), rStart,
0364                             kfOptions, tracks);
0365   if (!fitRes.ok()) {
0366     std::cout << "Fit failed" << std::endl;
0367     return ss.str();
0368   }
0369   auto& track = *fitRes;
0370
0371   // Draw the track
0372   std::cout << "Draw the fitted track" << std::endl;
0373   double momentumScale = 10;
0374   double localErrorScale = 100.;
0375   double directionErrorScale = 100000;
0376
0377   ViewConfig scolor{.color = {214, 214, 214}};
0378   ViewConfig mcolor{.color = {255, 145, 48}};
0379   mcolor.offset = -0.01;
0380   ViewConfig ppcolor{.color = {51, 204, 51}};
0381   ppcolor.offset = -0.02;
0382   ViewConfig fpcolor{.color = {255, 255, 0}};
0383   fpcolor.offset = -0.03;
0384   ViewConfig spcolor{.color = {0, 125, 255}};
0385   spcolor.offset = -0.04;
0386
0387   EventDataView3D::drawMultiTrajectory(
0388       helper, tracks.trackStateContainer(), track.tipIndex(), tgContext,
0389       momentumScale, localErrorScale, directionErrorScale, scolor, mcolor,
0390       ppcolor, fpcolor, spcolor);
0391
0392   helper.write("EventData_MultiTrajectory");
0393   helper.write(ss);
0394
0395   return ss.str();
0396 }
0397
0398 }  // namespace Acts::EventDataView3DTest