src/Root/MuonVisualization.cpp

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 "ActsExamples/Root/MuonVisualization.hpp"
0010
0011 #include "Acts/Definitions/Algebra.hpp"
0012 #include "Acts/Definitions/Units.hpp"
0013 #include "Acts/Geometry/GeometryIdentifier.hpp"
0014 #include "Acts/Geometry/TrackingVolume.hpp"
0015 #include "Acts/Surfaces/LineBounds.hpp"
0016 #include "Acts/Surfaces/RectangleBounds.hpp"
0017 #include "Acts/Surfaces/Surface.hpp"
0018 #include "Acts/Surfaces/SurfaceBounds.hpp"
0019 #include "Acts/Surfaces/TrapezoidBounds.hpp"
0020 #include "Acts/Utilities/ArrayHelpers.hpp"
0021 #include "Acts/Utilities/Logger.hpp"
0022 #include "Acts/Utilities/RangeXD.hpp"
0023 #include "Acts/Utilities/StringHelpers.hpp"
0024 #include "ActsExamples/EventData/SimParticle.hpp"
0025
0026 #include <algorithm>
0027 #include <cassert>
0028 #include <format>
0029 #include <memory>
0030 #include <mutex>
0031 #include <vector>
0032
0033 #include "TArrow.h"
0034 #include "TBox.h"
0035 #include "TCanvas.h"
0036 #include "TEllipse.h"
0037 #include "TH2D.h"
0038 #include "TH2I.h"
0039 #include "TLegend.h"
0040 #include "TMarker.h"
0041 #include "TObject.h"
0042 #include "TStyle.h"
0043
0044 using namespace Acts;
0045 using namespace Acts::UnitLiterals;
0046
0047 namespace {
0048
0049 /// @brief Returns the chamber ID from a full geometry identifier
0050 GeometryIdentifier toChamberId(const GeometryIdentifier& id) {
0051   return GeometryIdentifier{}.withVolume(id.volume()).withLayer(id.layer());
0052 }
0053
0054 /// @brief Returns the half-height of a trapezoid / rectangular bounds
0055 double halfHeight(const SurfaceBounds& bounds) {
0056   if (bounds.type() == SurfaceBounds::BoundsType::eRectangle) {
0057     return static_cast<const RectangleBounds&>(bounds).get(
0058         RectangleBounds::eMaxY);
0059   }
0060   // Trapezoid -> endcap
0061   else if (bounds.type() == SurfaceBounds::BoundsType::eTrapezoid) {
0062     return static_cast<const TrapezoidBounds&>(bounds).get(
0063         TrapezoidBounds::eHalfLengthY);
0064   }
0065   return std::numeric_limits<double>::max();
0066 }
0067
0068 /// @brief Draw a circle at position
0069 std::unique_ptr<TEllipse> drawCircle(const double x, double y, const double r,
0070                                      const int color = kBlack,
0071                                      const int fillStyle = 0) {
0072   auto circle = std::make_unique<TEllipse>(x, y, r);
0073   circle->SetLineColor(color);
0074   circle->SetFillStyle(fillStyle);
0075   circle->SetLineWidth(1);
0076   circle->SetFillColorAlpha(color, 0.2);
0077   return circle;
0078 }
0079
0080 /// @brief Draw a box at position
0081 std::unique_ptr<TBox> drawBox(const double cX, const double cY, const double wX,
0082                               const double wY, const int color = kBlack,
0083                               const int fillStyle = 3344) {
0084   auto box = std::make_unique<TBox>(cX - 0.5 * wX, cY - 0.5 * wY, cX + 0.5 * wX,
0085                                     cY + 0.5 * wY);
0086   box->SetLineColor(color);
0087   box->SetFillStyle(fillStyle);
0088   box->SetLineWidth(1);
0089   box->SetFillColorAlpha(color, 0.2);
0090   return box;
0091 }
0092
0093 /// @brief Determines the axis ranges in the z-y plane to draw
0094 ///        the measurements on the canvas
0095 RangeXD<2, double> canvasRanges(
0096     const ActsExamples::MuonSpacePointBucket& bucket) {
0097   RangeXD<2, double> ranges{
0098       filledArray<double, 2>(std::numeric_limits<double>::max()),
0099       filledArray<double, 2>(-std::numeric_limits<double>::max())};
0100   // Extra margin such that the canvas axes don't overlap with the depicted
0101   // measurements
0102   constexpr double extra = 3._cm;
0103   for (const auto& sp : bucket) {
0104     const Vector3& pos = sp.localPosition();
0105     ranges.expand(0, pos.z() - extra, pos.z() + extra);
0106     ranges.expand(1, pos.y() - extra, pos.y() + extra);
0107   }
0108   return ranges;
0109 }
0110
0111 /// @brief Returns whether a surface should be drawn on the canvas
0112 bool isSurfaceToDraw(const Surface& surface,
0113                      const RangeXD<2, double>& canvasBoundaries,
0114                      const Vector3& posInFrame) {
0115   // Draw only active surfaces
0116   if (surface.surfacePlacement() == nullptr) {
0117     return false;
0118   }
0119   const auto& bounds = surface.bounds();
0120
0121   if (surface.type() == Surface::SurfaceType::Plane) {
0122     const double hL{halfHeight(bounds)};
0123     // check whether the surface is inside the visible range
0124     const double minZ = std::max(posInFrame.z() - hL, canvasBoundaries.min(0));
0125     const double maxZ = std::min(posInFrame.z() + hL, canvasBoundaries.max(0));
0126     // The maximum is below the left side of the strip plane
0127     // or the minimum is above the right side of the plane
0128     return maxZ > posInFrame.z() - hL && minZ < posInFrame.z() + hL &&
0129            posInFrame.y() > canvasBoundaries.min(1) &&
0130            posInFrame.y() < canvasBoundaries.max(1);
0131   } else if (surface.type() == Surface::SurfaceType::Straw) {
0132     const double r = static_cast<const LineBounds&>(bounds).get(LineBounds::eR);
0133     // Check that the straw surface is well embedded on the canvas
0134     return posInFrame.y() - r > canvasBoundaries.min(1) &&
0135            posInFrame.y() + r < canvasBoundaries.max(1) &&
0136            posInFrame.z() - r > canvasBoundaries.min(0) &&
0137            posInFrame.z() + r < canvasBoundaries.max(0);
0138   }
0139
0140   return false;
0141 }
0142
0143 }  // namespace
0144
0145 namespace ActsExamples {
0146
0147 void visualizeMuonSpacePoints(const std::string& outputPath,
0148                               const GeometryContext& gctx,
0149                               const MuonSpacePointBucket& bucket,
0150                               const SimHitContainer& simHits,
0151                               const SimParticleContainer& simParticles,
0152                               const TrackingGeometry& trackingGeometry,
0153                               const Acts::Logger& logger) {
0154   // Helper function to transform from hit frame to space point frame
0155   auto toSpacePointFrame = [&](const GeometryIdentifier& hitId) -> Transform3 {
0156     const Surface* hitSurf = trackingGeometry.findSurface(hitId);
0157     assert(hitSurf != nullptr);
0158
0159     // Fetch the parent volume to express all points in the same coordinate
0160     // system
0161     const TrackingVolume* volume =
0162         trackingGeometry.findVolume(toChamberId(hitId));
0163     assert(volume != nullptr);
0164
0165     // Transformation to the common coordinate system of all space points
0166     const Transform3 parentTrf{AngleAxis3{90._degree, Vector3::UnitZ()} *
0167                                volume->globalToLocalTransform(gctx) *
0168                                hitSurf->localToGlobalTransform(gctx)};
0169     ACTS_VERBOSE("Transform into space point frame for surface "
0170                  << hitId << " is \n"
0171                  << Acts::toString(parentTrf));
0172     return parentTrf;
0173   };
0174
0175   auto canvas = std::make_unique<TCanvas>("can", "can", 600, 600);
0176   canvas->cd();
0177
0178   const GeometryIdentifier chambId = toChamberId(bucket.front().geometryId());
0179
0180   std::vector<std::unique_ptr<TObject>> primitives{};
0181
0182   const RangeXD<2, double> canvasBound{canvasRanges(bucket)};
0183   /// Draw the frame
0184   auto frame = std::make_unique<TH2I>("frame", "frame;z [mm];y [mm]", 1,
0185                                       canvasBound.min(0), canvasBound.max(0), 1,
0186                                       canvasBound.min(1), canvasBound.max(1));
0187   frame->Draw("AXIS");
0188
0189   // Loop over all surfaces inside the chamber volume to draw the ones covered
0190   // by the canvas
0191   const TrackingVolume* chambVolume = trackingGeometry.findVolume(chambId);
0192   assert(chambVolume != nullptr);
0193   chambVolume->apply(overloaded{[&](const Surface& surface) {
0194     const Vector3 pos = toSpacePointFrame(surface.geometryId()).translation();
0195     if (!isSurfaceToDraw(surface, canvasBound, pos)) {
0196       return;
0197     }
0198     const auto& bounds = surface.bounds();
0199     if (surface.type() == Surface::SurfaceType::Plane) {
0200       const double hL{halfHeight(bounds)};
0201       const double minZ = std::max(pos.z() - hL, canvasBound.min(0));
0202       const double maxZ = std::min(pos.z() + hL, canvasBound.max(0));
0203       primitives.push_back(drawBox(0.5 * (minZ + maxZ), pos.y(), maxZ - minZ,
0204                                    0.3_cm, kBlack, 0));
0205
0206     } else if (surface.type() == Surface::SurfaceType::Straw) {
0207       const double r =
0208           static_cast<const LineBounds&>(bounds).get(LineBounds::eR);
0209       primitives.push_back(drawCircle(pos.z(), pos.y(), r, kBlack, 0));
0210     }
0211   }});
0212
0213   for (auto& sp : bucket) {
0214     const Vector3& pos = sp.localPosition();
0215     if (sp.isStraw()) {
0216       primitives.push_back(
0217           drawCircle(pos.z(), pos.y(), sp.driftRadius(), kRed, 0));
0218     } else {
0219       primitives.push_back(
0220           drawBox(pos.z(), pos.y(), 3._cm, 0.5_cm, kRed, 1001));
0221     }
0222   }
0223   // Finally draw the muon trajectory
0224   for (const auto& simHit : simHits) {
0225     if (chambId != toChamberId(simHit.geometryId())) {
0226       continue;
0227     }
0228     const auto simPartItr = simParticles.find(simHit.particleId());
0229     if (simPartItr == simParticles.end() ||
0230         (*simPartItr).hypothesis() != ParticleHypothesis::muon()) {
0231       continue;
0232     }
0233     const auto toSpTrf = toSpacePointFrame(simHit.geometryId()) *
0234                          trackingGeometry.findSurface(simHit.geometryId())
0235                              ->localToGlobalTransform(gctx)
0236                              .inverse();
0237     const Vector3 pos = toSpTrf * simHit.position();
0238     const Vector3 dir = toSpTrf.linear() * simHit.direction();
0239     constexpr double arrowL = 1._cm;
0240     const Vector3 start = pos - 0.5 * arrowL * dir;
0241     const Vector3 end = pos + 0.5 * arrowL * dir;
0242     auto arrow =
0243         std::make_unique<TArrow>(start.z(), start.y(), end.z(), end.y(), 0.03);
0244     arrow->SetLineColor(kBlue + 1);
0245     arrow->SetLineWidth(1);
0246     arrow->SetLineStyle(kSolid);
0247     primitives.push_back(std::move(arrow));
0248   }
0249
0250   for (auto& prim : primitives) {
0251     prim->Draw();
0252   }
0253   canvas->SaveAs(outputPath.c_str());
0254 }
0255
0256 void visualizeMuonHoughMaxima(
0257     const std::string& outputPath, const MuonSpacePoint::MuonId& bucketId,
0258     const std::vector<Acts::HoughTransformUtils::PeakFinders::IslandsAroundMax<
0259         const MuonSpacePoint*>::Maximum>& maxima,
0260     const Acts::HoughTransformUtils::HoughPlane<const MuonSpacePoint*>& plane,
0261     const Acts::HoughTransformUtils::HoughAxisRanges& axis,
0262     const MuonSegmentContainer& truthSegments, const Acts::Logger& logger) {
0263   static std::mutex canvasMutex{};
0264   std::lock_guard guard{canvasMutex};
0265
0266   TCanvas canvas("houghCanvas", "", 800, 800);
0267   canvas.SetRightMargin(0.12);
0268   canvas.SetLeftMargin(0.12);
0269   gStyle->SetPalette(kGreyScale);
0270   gStyle->SetOptStat(0);
0271
0272   std::vector<std::unique_ptr<TObject>> primitives;
0273
0274   /// Save the hough accumulator as histogram
0275   TH2D houghHistoForPlot("houghHist", "HoughPlane;tan(#alpha);z0 [mm]",
0276                          plane.nBinsX(), axis.xMin, axis.xMax, plane.nBinsY(),
0277                          axis.yMin, axis.yMax);
0278   houghHistoForPlot.SetTitle(
0279       std::format("Station {:}, side {:}, sector {:2d}",
0280                   MuonSpacePoint::MuonId::toString(bucketId.msStation()),
0281                   MuonSpacePoint::MuonId::toString(bucketId.side()),
0282                   bucketId.sector())
0283           .c_str());
0284
0285   /** Copy the plane content into the histogram */
0286   for (int bx = 0; bx < houghHistoForPlot.GetNbinsX(); ++bx) {
0287     for (int by = 0; by < houghHistoForPlot.GetNbinsY(); ++by) {
0288       houghHistoForPlot.SetBinContent(bx + 1, by + 1, plane.nHits(bx, by));
0289     }
0290   }
0291   /** Set the contours */
0292   auto maxHitsAsInt = static_cast<int>(plane.maxHits());
0293   houghHistoForPlot.SetContour(maxHitsAsInt + 1);
0294   for (int k = 0; k < maxHitsAsInt + 1; ++k) {
0295     houghHistoForPlot.SetContourLevel(k, k - 0.5);
0296   }
0297
0298   auto legend = std::make_unique<TLegend>(0.5, 0.7, 1. - gPad->GetRightMargin(),
0299                                           1. - gPad->GetTopMargin());
0300   legend->SetBorderSize(0);
0301   legend->SetFillStyle(0);
0302
0303   /** Fetch the true parameters */
0304   MuonSegmentContainer::const_iterator truthItr = truthSegments.begin();
0305   const int trueCoord = bucketId.measuresEta() ? Acts::eY : Acts::eX;
0306   while ((truthItr = std::find_if(truthItr, truthSegments.end(),
0307                                   [bucketId](const MuonSegment& seg) {
0308                                     return seg.id().sameStation(bucketId);
0309                                   })) != truthSegments.end()) {
0310     const MuonSegment& truthSeg{*truthItr};
0311     const float tanAlpha =
0312         truthSeg.localDirection()[trueCoord] / truthSeg.localDirection().z();
0313     const float intercept = truthSeg.localPosition()[trueCoord];
0314
0315     auto trueMarker =
0316         std::make_unique<TMarker>(tanAlpha, intercept, kOpenCrossX);
0317     trueMarker->SetMarkerSize(3);
0318     trueMarker->SetMarkerColor(kRed);
0319     legend->AddEntry(trueMarker.get(), "True coordinates");
0320     primitives.push_back(std::move(trueMarker));
0321     ACTS_VERBOSE("Draw true segment " << truthSeg.id()
0322                                       << " in bucket: " << bucketId);
0323     ++truthItr;
0324   }
0325
0326   bool addedLeg{false};
0327   for (const auto& max : maxima) {
0328     auto marker = std::make_unique<TMarker>(max.x, max.y, kFullSquare);
0329     marker->SetMarkerSize(1);
0330     marker->SetMarkerColor(kBlue);
0331
0332     auto box = std::make_unique<TBox>(max.x - max.wx, max.y - max.wy,
0333                                       max.x + max.wx, max.y + max.wy);
0334     box->SetLineColor(kBlue);
0335     box->SetFillStyle(1001);
0336     box->SetFillColorAlpha(kBlue, 0.1);
0337     box->SetLineWidth(0);
0338     if (!addedLeg) {
0339       legend->AddEntry(marker.get(), "Hough maxima");
0340       legend->AddEntry(box.get(), "Hough uncertainties");
0341       addedLeg = true;
0342     }
0343     primitives.push_back(std::move(box));
0344     primitives.push_back(std::move(marker));
0345   }
0346   primitives.emplace_back(std::move(legend));
0347   for (auto& prim : primitives) {
0348     prim->Draw();
0349   }
0350
0351   canvas.SaveAs(outputPath.c_str());
0352 }
0353
0354 }  // namespace ActsExamples