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