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