Back to home page

sPhenix code displayed by LXR

 
 

    


File indexing completed on 2026-07-16 08:08:01

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/Digitization/MuonSpacePointDigitizer.hpp"
0010 
0011 #include "Acts/Definitions/Units.hpp"
0012 #include "Acts/Geometry/GeometryContext.hpp"
0013 #include "Acts/Geometry/VolumeBounds.hpp"
0014 #include "Acts/Surfaces/LineBounds.hpp"
0015 #include "Acts/Surfaces/RectangleBounds.hpp"
0016 #include "Acts/Surfaces/TrapezoidBounds.hpp"
0017 #include "Acts/Surfaces/detail/LineHelper.hpp"
0018 #include "Acts/Surfaces/detail/PlanarHelper.hpp"
0019 #include "Acts/Utilities/ArrayHelpers.hpp"
0020 #include "Acts/Utilities/Helpers.hpp"
0021 #include "Acts/Utilities/MathHelpers.hpp"
0022 #include "Acts/Utilities/StringHelpers.hpp"
0023 #include "ActsExamples/Digitization/Smearers.hpp"
0024 #include "ActsExamples/EventData/MuonSpacePoint.hpp"
0025 
0026 #include <algorithm>
0027 #include <format>
0028 #include <map>
0029 #include <ranges>
0030 
0031 #include "TArrow.h"
0032 #include "TBox.h"
0033 #include "TCanvas.h"
0034 #include "TEllipse.h"
0035 #include "TH2I.h"
0036 #include "TROOT.h"
0037 #include "TStyle.h"
0038 
0039 using namespace Acts;
0040 using namespace detail::LineHelper;
0041 using namespace PlanarHelper;
0042 using namespace UnitLiterals;
0043 
0044 namespace {
0045 
0046 /// @brief Quanitze the hit position to a strip
0047 /// @param x: Actual traversing of the particle
0048 /// @param pitch: Pitch between two strips
0049 constexpr double quantize(const double x, const double pitch) {
0050   if (x >= 0.) {
0051     return std::max(std::floor(x - 0.5 * pitch) / pitch, 0.) * pitch;
0052   }
0053   return -quantize(-x, pitch);
0054 }
0055 /// @brief Returns the half-height of a trapezoid / rectangular bounds
0056 /// @param bounds: Rectangle / Trapezoid bounds to fetch the half height from
0057 double halfHeight(const SurfaceBounds& bounds) {
0058   if (bounds.type() == SurfaceBounds::BoundsType::eRectangle) {
0059     return static_cast<const RectangleBounds&>(bounds).get(
0060         RectangleBounds::eMaxY);
0061   }
0062   // Trapezoid -> endcap
0063   else if (bounds.type() == SurfaceBounds::BoundsType::eTrapezoid) {
0064     return static_cast<const TrapezoidBounds&>(bounds).get(
0065         TrapezoidBounds::eHalfLengthY);
0066   }
0067   return std::numeric_limits<double>::max();
0068 }
0069 /// @brief Draw a circle at position
0070 std::unique_ptr<TEllipse> drawCircle(const double x, double y, const double r,
0071                                      const int color = kBlack,
0072                                      const int fillStyle = 0) {
0073   auto circle = std::make_unique<TEllipse>(x, y, r);
0074   circle->SetLineColor(color);
0075   circle->SetFillStyle(fillStyle);
0076   circle->SetLineWidth(1);
0077   circle->SetFillColorAlpha(color, 0.2);
0078   return circle;
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 constexpr GeometryIdentifier toChamberId(const GeometryIdentifier& id) {
0094   return GeometryIdentifier{}.withVolume(id.volume()).withLayer(id.layer());
0095 }
0096 
0097 }  // namespace
0098 
0099 namespace ActsExamples {
0100 
0101 MuonSpacePointDigitizer::MuonSpacePointDigitizer(const Config& cfg,
0102                                                  Logging::Level lvl)
0103     : IAlgorithm("MuonSpacePointDigitizer", lvl), m_cfg{cfg} {
0104   if (m_cfg.inputSimHits.empty()) {
0105     throw std::invalid_argument("No sim hits have been parsed ");
0106   }
0107   if (m_cfg.inputParticles.empty()) {
0108     throw std::invalid_argument("No simulated particles were parsed");
0109   }
0110   if (m_cfg.outputSpacePoints.empty()) {
0111     throw std::invalid_argument("No output space points were defined");
0112   }
0113   ACTS_DEBUG("Retrieve sim hits and particles from "
0114              << m_cfg.inputSimHits << " & " << m_cfg.inputParticles);
0115   ACTS_DEBUG("Write produced space points to " << m_cfg.outputSpacePoints);
0116   m_inputSimHits.initialize(m_cfg.inputSimHits);
0117   m_inputParticles.initialize(m_cfg.inputParticles);
0118   m_outputSpacePoints.initialize(m_cfg.outputSpacePoints);
0119 }
0120 
0121 ProcessCode MuonSpacePointDigitizer::initialize() {
0122   using enum ActsExamples::ProcessCode;
0123   if (!m_cfg.trackingGeometry) {
0124     ACTS_ERROR("No tracking geometry was parsed");
0125     return ABORT;
0126   }
0127   if (!m_cfg.randomNumbers) {
0128     ACTS_ERROR("No random number generator was parsed");
0129     return ABORT;
0130   }
0131   MuonSpacePointCalibrator::Config calibCfg{};
0132   m_cfg.calibrator =
0133       std::make_unique<MuonSpacePointCalibrator>(calibCfg, logger().clone());
0134 
0135   gROOT->SetStyle("ATLAS");
0136   return SUCCESS;
0137 }
0138 
0139 const MuonSpacePointCalibrator& MuonSpacePointDigitizer::calibrator() const {
0140   assert(m_cfg.calibrator != nullptr);
0141   return *m_cfg.calibrator;
0142 }
0143 const TrackingGeometry& MuonSpacePointDigitizer::trackingGeometry() const {
0144   assert(m_cfg.trackingGeometry != nullptr);
0145   return *m_cfg.trackingGeometry;
0146 }
0147 Transform3 MuonSpacePointDigitizer::toSpacePointFrame(
0148     const GeometryContext& gctx, const GeometryIdentifier& hitId) const {
0149   const Surface* hitSurf = trackingGeometry().findSurface(hitId);
0150   assert(hitSurf != nullptr);
0151 
0152   /// Fetch the parent volume to express all points in the same coordinate
0153   /// system
0154   const TrackingVolume* volume =
0155       trackingGeometry().findVolume(toChamberId(hitId));
0156   assert(volume != nullptr);
0157   /// Transformation to the common coordinate system of all space points
0158   const Transform3 parentTrf{AngleAxis3{90._degree, Vector3::UnitZ()} *
0159                              volume->globalToLocalTransform(gctx) *
0160                              hitSurf->localToGlobalTransform(gctx)};
0161   ACTS_VERBOSE("Transform into space point frame for surface "
0162                << hitId << " is \n"
0163                << toString(parentTrf));
0164   return parentTrf;
0165 }
0166 ProcessCode MuonSpacePointDigitizer::execute(
0167     const AlgorithmContext& ctx) const {
0168   const SimHitContainer& gotSimHits = m_inputSimHits(ctx);
0169   const SimParticleContainer& simParticles = m_inputParticles(ctx);
0170   ACTS_DEBUG("Retrieved " << gotSimHits.size() << " hits & "
0171                           << simParticles.size() << " associated particles.");
0172 
0173   MuonSpacePointContainer outSpacePoints{};
0174 
0175   const auto gctx = Acts::GeometryContext::dangerouslyDefaultConstruct();
0176 
0177   using MuonId_t = MuonSpacePoint::MuonId;
0178   auto rndEngine = m_cfg.randomNumbers->spawnGenerator(ctx);
0179   /// temporary output container to group the hits per chamber volume
0180   std::map<GeometryIdentifier, MuonSpacePointBucket> spacePointsPerChamber{};
0181   std::unordered_map<GeometryIdentifier, double> strawTimes{};
0182   std::multimap<GeometryIdentifier, std::array<double, 3>> stripTimes{};
0183 
0184   ACTS_DEBUG("Starting loop over modules ...");
0185   for (const auto& simHitsGroup : groupByModule(gotSimHits)) {
0186     // Manual pair unpacking instead of using
0187     //   auto [moduleGeoId, moduleSimHits] : ...
0188     // otherwise clang on macos complains that it is unable to capture the local
0189     // binding in the lambda used for visiting the smearer below.
0190     Acts::GeometryIdentifier moduleGeoId = simHitsGroup.first;
0191     const auto& moduleSimHits = simHitsGroup.second;
0192 
0193     const Surface* hitSurf = trackingGeometry().findSurface(moduleGeoId);
0194     assert(hitSurf != nullptr);
0195 
0196     const Transform3& surfLocToGlob{hitSurf->localToGlobalTransform(gctx)};
0197 
0198     /// Transformation to the common coordinate system of all space points
0199     const Transform3 parentTrf{toSpacePointFrame(gctx, moduleGeoId)};
0200     /// Retrieve the bounds
0201     const auto& bounds = hitSurf->bounds();
0202 
0203     // Iterate over all simHits in a single module
0204     for (auto h = moduleSimHits.begin(); h != moduleSimHits.end(); ++h) {
0205       const auto& simHit = *h;
0206 
0207       // Convert the hit trajectory into local coordinates
0208       const Vector3 locPos = surfLocToGlob.inverse() * simHit.position();
0209       const Vector3 locDir =
0210           surfLocToGlob.inverse().linear() * simHit.direction();
0211 
0212       ACTS_DEBUG("Process hit: " << toString(locPos)
0213                                  << ", dir: " << toString(locDir)
0214                                  << " recorded in a " << hitSurf->type()
0215                                  << " surface with id: " << moduleGeoId
0216                                  << ", bounds: " << bounds);
0217       bool convertSp{true};
0218 
0219       MuonSpacePoint newSp{};
0220       newSp.setGeometryId(moduleGeoId);
0221 
0222       const auto& calibCfg = calibrator().config();
0223       switch (hitSurf->type()) {
0224         /// Strip measurements
0225         using enum Surface::SurfaceType;
0226         case Plane: {
0227           ACTS_VERBOSE("Hit is recorded in a strip detector ");
0228           auto planeCross =
0229               intersectPlane(locPos, locDir, Vector3::UnitZ(), 0.);
0230           const auto hitPos = planeCross.position();
0231           Vector3 smearedHit{Vector3::Zero()};
0232           switch (bounds.type()) {
0233             case SurfaceBounds::BoundsType::eRectangle: {
0234               smearedHit[ePos0] =
0235                   quantize(hitPos[ePos0], calibCfg.rpcPhiStripPitch);
0236               smearedHit[ePos1] =
0237                   quantize(hitPos[ePos1], calibCfg.rpcEtaStripPitch);
0238               ACTS_VERBOSE("Position before "
0239                            << toString(hitPos) << ", after smearing "
0240                            << toString(smearedHit) << ", " << bounds);
0241 
0242               if (!bounds.inside(
0243                       Vector2{smearedHit[ePos0], smearedHit[ePos1]})) {
0244                 convertSp = false;
0245                 break;
0246               }
0247               auto ranges = stripTimes.equal_range(moduleGeoId);
0248               for (auto digitHitItr = ranges.first;
0249                    digitHitItr != ranges.second; ++digitHitItr) {
0250                 const auto& existCoords = digitHitItr->second;
0251                 /// Same virtual strip point is digitized
0252                 if (std::abs(existCoords[0] - smearedHit[ePos0]) <
0253                         Acts::s_epsilon &&
0254                     std::abs(existCoords[1] - smearedHit[ePos1]) <
0255                         Acts::s_epsilon &&
0256                     simHit.time() - existCoords[2] < config().rpcDeadTime) {
0257                   convertSp = false;
0258                   break;
0259                 }
0260                 if (!convertSp) {
0261                   break;
0262                 }
0263               }
0264               /// Mark that a new hit has been recorded at this position & time
0265               /// Subsequent hits are rejected if they remain within the dead
0266               /// time
0267               stripTimes.insert(std::make_pair(
0268                   moduleGeoId, std::array{smearedHit[ePos0], smearedHit[ePos1],
0269                                           simHit.time()}));
0270 
0271               /// Time digitization
0272               if (config().digitizeTime) {
0273                 assert(calibCfg.rpcTimeResolution > 0.);
0274                 const double stripTime =
0275                     (*Digitization::Gauss{calibCfg.rpcTimeResolution}(
0276                          simHit.time(), rndEngine))
0277                         .first;
0278                 newSp.setTime(stripTime);
0279               }
0280               newSp.setCovariance(
0281                   calibCfg.rpcPhiStripPitch, calibCfg.rpcEtaStripPitch,
0282                   m_cfg.digitizeTime ? calibCfg.rpcTimeResolution : 0.);
0283 
0284               break;
0285             }
0286             /// Endcap strips not yet available
0287             case SurfaceBounds::BoundsType::eTrapezoid:
0288               break;
0289             default:
0290               convertSp = false;
0291           }
0292           /// Define the space point coordinates
0293           if (convertSp) {
0294             newSp.defineCoordinates(
0295                 Vector3{parentTrf * smearedHit},
0296                 Vector3{parentTrf.linear().col(Acts::ePos1)},
0297                 Vector3{parentTrf.linear().col(Acts::ePos0)});
0298             MuonId_t id{};
0299             /// @todo Refine me using the volume name
0300             id.setChamber(MuonId_t::StationName::BIS,
0301                           simHit.position().z() > 0 ? MuonId_t::DetSide::A
0302                                                     : MuonId_t::DetSide::C,
0303                           1, MuonId_t::TechField::Rpc);
0304             id.setCoordFlags(true, true);
0305             newSp.setId(id);
0306           }
0307 
0308           break;
0309         }
0310         case Straw: {
0311           auto closeApproach = lineIntersect<3>(
0312               Vector3::Zero(), Vector3::UnitZ(), locPos, locDir);
0313           const auto nominalPos = closeApproach.position();
0314           const double unsmearedR = fastHypot(nominalPos.x(), nominalPos.y());
0315           ACTS_VERBOSE("Hit is recorded in a straw detector, R: "
0316                        << unsmearedR << ", " << bounds);
0317 
0318           const double uncert = calibrator().driftRadiusUncert(unsmearedR);
0319           /// Reject unsmearable hits
0320           if (uncert <= std::numeric_limits<double>::epsilon()) {
0321             convertSp = false;
0322             break;
0323           }
0324           double driftR =
0325               (*Digitization::Gauss{uncert}(unsmearedR, rndEngine)).first;
0326 
0327           // bounds
0328           const auto& lBounds = static_cast<const LineBounds&>(bounds);
0329           const double maxR = lBounds.get(LineBounds::eR);
0330           const double maxZ = lBounds.get(LineBounds::eHalfLengthZ);
0331           /// The generated hit is unphysical
0332           if (driftR < 0. || driftR > maxR || std::abs(nominalPos.z()) > maxZ) {
0333             convertSp = false;
0334             break;
0335           }
0336           if (auto insertItr =
0337                   strawTimes.insert(std::make_pair(moduleGeoId, simHit.time()));
0338               !insertItr.second) {
0339             if (simHit.time() - insertItr.first->second > m_cfg.strawDeadTime) {
0340               insertItr.first->second = simHit.time();
0341             } else {
0342               convertSp = false;
0343               break;
0344             }
0345           }
0346 
0347           const double sigmaZ = 0.5 * maxZ;
0348 
0349           newSp.setRadius(driftR);
0350           newSp.setCovariance(square(uncert), square(sigmaZ), 0.);
0351 
0352           newSp.defineCoordinates(
0353               Vector3{parentTrf.translation()},
0354               Vector3{parentTrf.linear() * Vector3::UnitZ()},
0355               Vector3{parentTrf.linear() * Vector3::UnitX()});
0356           MuonId_t id{};
0357           /// @todo Refine me using the volume name
0358           id.setChamber(MuonId_t::StationName::BIS,
0359                         simHit.position().z() > 0 ? MuonId_t::DetSide::A
0360                                                   : MuonId_t::DetSide::C,
0361                         1, MuonId_t::TechField::Mdt);
0362           id.setCoordFlags(true, false);
0363           newSp.setId(id);
0364 
0365           break;
0366         }
0367         default:
0368           ACTS_DEBUG(
0369               "Unsupported detector case in muon space point digitizer.");
0370           convertSp = false;
0371       }
0372 
0373       if (convertSp) {
0374         spacePointsPerChamber[toChamberId(moduleGeoId)].push_back(
0375             std::move(newSp));
0376       }
0377     }
0378 
0379     for (auto& [volId, bucket] : spacePointsPerChamber) {
0380       std::ranges::sort(bucket,
0381                         [](const MuonSpacePoint& a, const MuonSpacePoint& b) {
0382                           return a.localPosition().z() < b.localPosition().z();
0383                         });
0384       if (logger().doPrint(Logging::Level::VERBOSE)) {
0385         std::stringstream sstr{};
0386         for (const auto& spacePoint : bucket) {
0387           sstr << " *** " << spacePoint << std::endl;
0388         }
0389         ACTS_VERBOSE("Safe " << bucket.size() << " space points for chamber "
0390                              << volId << "\n"
0391                              << sstr.str());
0392       }
0393       visualizeBucket(ctx, gctx, bucket);
0394       outSpacePoints.push_back(std::move(bucket));
0395     }
0396   }
0397 
0398   m_outputSpacePoints(ctx, std::move(outSpacePoints));
0399 
0400   return ProcessCode::SUCCESS;
0401 }
0402 
0403 RangeXD<2, double> MuonSpacePointDigitizer::canvasRanges(
0404     const MuonSpacePointBucket& bucket) const {
0405   RangeXD<2, double> ranges{
0406       filledArray<double, 2>(std::numeric_limits<double>::max()),
0407       filledArray<double, 2>(-std::numeric_limits<double>::max())};
0408   // Extra margin such that the canvas axes don't overlap with the depicted
0409   // measurements
0410   constexpr double extra = 3._cm;
0411   for (const auto& sp : bucket) {
0412     const Vector3& pos = sp.localPosition();
0413     ranges.expand(0, pos.z() - extra, pos.z() + extra);
0414     ranges.expand(1, pos.y() - extra, pos.y() + extra);
0415   }
0416   return ranges;
0417 }
0418 
0419 bool MuonSpacePointDigitizer::isSurfaceToDraw(
0420     const Acts::GeometryContext& gctx, const Surface& surface,
0421     const RangeXD<2, double>& canvasBoundaries) const {
0422   // Draw only active surfaces
0423   if (!surface.isSensitive()) {
0424     return false;
0425   }
0426   // surface position in the frame
0427   const Vector3 pos =
0428       toSpacePointFrame(gctx, surface.geometryId()).translation();
0429   const auto& bounds = surface.bounds();
0430 
0431   if (surface.type() == Surface::SurfaceType::Plane) {
0432     const double hL{halfHeight(bounds)};
0433     // check whether the surface is inside the visible range
0434     const double minZ = std::max(pos.z() - hL, canvasBoundaries.min(0));
0435     const double maxZ = std::min(pos.z() + hL, canvasBoundaries.max(0));
0436     // The maximum is below the left side of the strip plane
0437     // or the minimum is above the right side of the plane
0438     return maxZ > pos.z() - hL && minZ < pos.z() + hL &&
0439            pos.y() > canvasBoundaries.min(1) &&
0440            pos.y() < canvasBoundaries.max(1);
0441   } else if (surface.type() == Surface::SurfaceType::Straw) {
0442     const double r = static_cast<const LineBounds&>(bounds).get(LineBounds::eR);
0443     // Check that the straw surface is well embedded on the canvas
0444     return pos.y() - r > canvasBoundaries.min(1) &&
0445            pos.y() + r < canvasBoundaries.max(1) &&
0446            pos.z() - r > canvasBoundaries.min(0) &&
0447            pos.z() + r < canvasBoundaries.max(0);
0448   }
0449 
0450   return false;
0451 }
0452 void MuonSpacePointDigitizer::visualizeBucket(
0453     const AlgorithmContext& ctx, const GeometryContext& gctx,
0454     const MuonSpacePointBucket& bucket) const {
0455   if (!m_cfg.dumpVisualization) {
0456     return;
0457   }
0458   auto canvas = std::make_unique<TCanvas>("can", "can", 600, 600);
0459   canvas->cd();
0460 
0461   const GeometryIdentifier chambId = toChamberId(bucket.front().geometryId());
0462 
0463   std::vector<std::unique_ptr<TObject>> primitives{};
0464 
0465   const RangeXD<2, double> canvasBound{canvasRanges(bucket)};
0466   /// Draw the frame
0467   auto frame = std::make_unique<TH2I>("frame", "frame;z [mm];y [mm]", 1,
0468                                       canvasBound.min(0), canvasBound.max(0), 1,
0469                                       canvasBound.min(1), canvasBound.max(1));
0470   frame->Draw("AXIS");
0471 
0472   // Loop over all surfaces inside the chamber volume to draw the ones covered
0473   // by the canvas
0474   const TrackingVolume* chambVolume = trackingGeometry().findVolume(chambId);
0475   assert(chambVolume != nullptr);
0476   chambVolume->apply(overloaded{
0477       [this, &canvasBound, &gctx, &primitives](const Surface& surface) {
0478         if (!isSurfaceToDraw(gctx, surface, canvasBound)) {
0479           return;
0480         }
0481         const Vector3 pos =
0482             toSpacePointFrame(gctx, surface.geometryId()).translation();
0483         const auto& bounds = surface.bounds();
0484         if (surface.type() == Surface::SurfaceType::Plane) {
0485           const double hL{halfHeight(bounds)};
0486           const double minZ = std::max(pos.z() - hL, canvasBound.min(0));
0487           const double maxZ = std::min(pos.z() + hL, canvasBound.max(0));
0488           primitives.push_back(drawBox(0.5 * (minZ + maxZ), pos.y(),
0489                                        maxZ - minZ, 0.3_cm, kBlack, 0));
0490 
0491         } else if (surface.type() == Surface::SurfaceType::Straw) {
0492           const double r =
0493               static_cast<const LineBounds&>(bounds).get(LineBounds::eR);
0494           primitives.push_back(drawCircle(pos.z(), pos.y(), r, kBlack, 0));
0495         }
0496       }});
0497 
0498   for (auto& sp : bucket) {
0499     const Vector3& pos = sp.localPosition();
0500     if (sp.isStraw()) {
0501       primitives.push_back(
0502           drawCircle(pos.z(), pos.y(), sp.driftRadius(), kRed, 0));
0503     } else {
0504       primitives.push_back(
0505           drawBox(pos.z(), pos.y(), 3._cm, 0.5_cm, kRed, 1001));
0506     }
0507   }
0508   // Finally draw the muon trajectory
0509   const SimHitContainer& gotSimHits = m_inputSimHits(ctx);
0510   const SimParticleContainer& simParticles = m_inputParticles(ctx);
0511   for (const auto& simHit : gotSimHits) {
0512     if (chambId != toChamberId(simHit.geometryId())) {
0513       continue;
0514     }
0515     const auto simPartItr = simParticles.find(simHit.particleId());
0516     if (simPartItr == simParticles.end() ||
0517         (*simPartItr).hypothesis() != ParticleHypothesis::muon()) {
0518       continue;
0519     }
0520     const auto toSpTrf = toSpacePointFrame(gctx, simHit.geometryId()) *
0521                          trackingGeometry()
0522                              .findSurface(simHit.geometryId())
0523                              ->localToGlobalTransform(gctx)
0524                              .inverse();
0525     const Vector3 pos = toSpTrf * simHit.position();
0526     const Vector3 dir = toSpTrf.linear() * simHit.direction();
0527     constexpr double arrowL = 1._cm;
0528     const Vector3 start = pos - 0.5 * arrowL * dir;
0529     const Vector3 end = pos + 0.5 * arrowL * dir;
0530     auto arrow =
0531         std::make_unique<TArrow>(start.z(), start.y(), end.z(), end.y(), 0.03);
0532     arrow->SetLineColor(kBlue + 1);
0533     arrow->SetLineWidth(1);
0534     arrow->SetLineStyle(kSolid);
0535     primitives.push_back(std::move(arrow));
0536   }
0537 
0538   for (auto& prim : primitives) {
0539     prim->Draw();
0540   }
0541   canvas->SaveAs(
0542       std::format("Event_{}_{}.pdf", ctx.eventNumber, chambVolume->volumeName())
0543           .c_str());
0544 }
0545 
0546 }  // namespace ActsExamples