Digitization/src/MuonSpacePointDigitizer.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/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
0105 ProcessCode MuonSpacePointDigitizer::initialize() {
0106   using enum ActsExamples::ProcessCode;
0107   if (!m_cfg.trackingGeometry) {
0108     ACTS_ERROR("No tracking geometry was parsed");
0109     return ABORT;
0110   }
0111   if (!m_cfg.randomNumbers) {
0112     ACTS_ERROR("No random number generator was parsed");
0113     return ABORT;
0114   }
0115   MuonSpacePointCalibrator::Config calibCfg{};
0116   m_cfg.calibrator =
0117       std::make_unique<MuonSpacePointCalibrator>(calibCfg, logger().clone());
0118
0119   if (m_cfg.inputSimHits.empty()) {
0120     ACTS_ERROR("No sim hits have been parsed ");
0121     return ABORT;
0122   }
0123   if (m_cfg.inputParticles.empty()) {
0124     ACTS_ERROR("No simulated particles were parsed");
0125     return ABORT;
0126   }
0127   if (m_cfg.outputSpacePoints.empty()) {
0128     ACTS_ERROR("No output space points were defined");
0129     return ABORT;
0130   }
0131   ACTS_DEBUG("Retrieve sim hits and particles from "
0132              << m_cfg.inputSimHits << " & " << m_cfg.inputParticles);
0133   ACTS_DEBUG("Write produced space points to " << m_cfg.outputSpacePoints);
0134   m_inputSimHits.initialize(m_cfg.inputSimHits);
0135   m_inputParticles.initialize(m_cfg.inputParticles);
0136   m_outputSpacePoints.initialize(m_cfg.outputSpacePoints);
0137   gROOT->SetStyle("ATLAS");
0138   return SUCCESS;
0139 }
0140
0141 const MuonSpacePointCalibrator& MuonSpacePointDigitizer::calibrator() const {
0142   assert(m_cfg.calibrator != nullptr);
0143   return *m_cfg.calibrator;
0144 }
0145 const TrackingGeometry& MuonSpacePointDigitizer::trackingGeometry() const {
0146   assert(m_cfg.trackingGeometry != nullptr);
0147   return *m_cfg.trackingGeometry;
0148 }
0149 Transform3 MuonSpacePointDigitizer::toSpacePointFrame(
0150     const GeometryContext& gctx, const GeometryIdentifier& hitId) const {
0151   const Surface* hitSurf = trackingGeometry().findSurface(hitId);
0152   assert(hitSurf != nullptr);
0153
0154   /// Fetch the parent volume to express all points in the same coordinate
0155   /// system
0156   const TrackingVolume* volume =
0157       trackingGeometry().findVolume(toChamberId(hitId));
0158   assert(volume != nullptr);
0159   /// Transformation to the common coordinate system of all space points
0160   const Transform3 parentTrf{AngleAxis3{90._degree, Vector3::UnitZ()} *
0161                              volume->itransform() * hitSurf->transform(gctx)};
0162   ACTS_VERBOSE("Transform into space point frame for surface "
0163                << hitId << " is \n"
0164                << toString(parentTrf));
0165   return parentTrf;
0166 }
0167 ProcessCode MuonSpacePointDigitizer::execute(
0168     const AlgorithmContext& ctx) const {
0169   const SimHitContainer& gotSimHits = m_inputSimHits(ctx);
0170   const SimParticleContainer& simParticles = m_inputParticles(ctx);
0171   ACTS_DEBUG("Retrieved " << gotSimHits.size() << " hits & "
0172                           << simParticles.size() << " associated particles.");
0173
0174   MuonSpacePointContainer outSpacePoints{};
0175
0176   GeometryContext gctx{};
0177
0178   using MuonId_t = MuonSpacePoint::MuonId;
0179   auto rndEngine = m_cfg.randomNumbers->spawnGenerator(ctx);
0180   /// temporary output container to group the hits per chamber volume
0181   std::map<GeometryIdentifier, MuonSpacePointBucket> spacePointsPerChamber{};
0182   std::unordered_map<GeometryIdentifier, double> strawTimes{};
0183   std::multimap<GeometryIdentifier, std::array<double, 3>> stripTimes{};
0184
0185   for (const auto& hit : gotSimHits) {
0186     const GeometryIdentifier hitId = hit.geometryId();
0187
0188     const Surface* hitSurf = trackingGeometry().findSurface(hitId);
0189     assert(hitSurf != nullptr);
0190
0191     const Transform3& surfLocToGlob{hitSurf->transform(gctx)};
0192
0193     // Convert the hit trajectory into local coordinates
0194     const Vector3 locPos = surfLocToGlob.inverse() * hit.position();
0195     const Vector3 locDir = surfLocToGlob.inverse().linear() * hit.direction();
0196
0197     const auto& bounds = hitSurf->bounds();
0198     ACTS_DEBUG("Process hit: "
0199                << toString(locPos) << ", dir: " << toString(locDir)
0200                << "recorded in a "
0201                << Surface::s_surfaceTypeNames[toUnderlying(hitSurf->type())]
0202                << " surface with id: " << hitId << ", bounds: " << bounds);
0203     bool convertSp{true};
0204
0205     MuonSpacePoint newSp{};
0206     newSp.setGeometryId(hitId);
0207
0208     /// Transformation to the common coordinate system of all space points
0209     const Transform3 parentTrf{toSpacePointFrame(gctx, hitId)};
0210
0211     const auto& calibCfg = calibrator().config();
0212     switch (hitSurf->type()) {
0213       /// Strip measurements
0214       using enum Surface::SurfaceType;
0215       case Plane: {
0216         ACTS_VERBOSE("Hit is recorded in a strip detector ");
0217         auto planeCross = intersectPlane(locPos, locDir, Vector3::UnitZ(), 0.);
0218         const auto hitPos = planeCross.position();
0219         Vector3 smearedHit{Vector3::Zero()};
0220         switch (bounds.type()) {
0221           case SurfaceBounds::BoundsType::eRectangle: {
0222             smearedHit[ePos0] =
0223                 quantize(hitPos[ePos0], calibCfg.rpcPhiStripPitch);
0224             smearedHit[ePos1] =
0225                 quantize(hitPos[ePos1], calibCfg.rpcEtaStripPitch);
0226             ACTS_VERBOSE("Position before "
0227                          << toString(hitPos) << ", after smearing"
0228                          << toString(smearedHit) << ", " << bounds);
0229
0230             if (!bounds.inside(Vector2{smearedHit[ePos0], smearedHit[ePos1]})) {
0231               convertSp = false;
0232               break;
0233             }
0234             auto ranges = stripTimes.equal_range(hitId);
0235             for (auto digitHitItr = ranges.first; digitHitItr != ranges.second;
0236                  ++digitHitItr) {
0237               const auto& existCoords = digitHitItr->second;
0238               /// Same virtual strip point is digitized
0239               if (std::abs(existCoords[0] - smearedHit[ePos0]) <
0240                       std::numeric_limits<double>::epsilon() &&
0241                   std::abs(existCoords[1] - smearedHit[ePos1]) <
0242                       std::numeric_limits<double>::epsilon() &&
0243                   hit.time() - existCoords[2] < config().rpcDeadTime) {
0244                 convertSp = false;
0245                 break;
0246               }
0247               if (!convertSp) {
0248                 break;
0249               }
0250             }
0251             /// Mark the
0252             stripTimes.insert(std::make_pair(
0253                 hitId,
0254                 std::array{smearedHit[ePos0], smearedHit[ePos1], hit.time()}));
0255
0256             /// Time digitization
0257             if (config().digitizeTime) {
0258               assert(calibCfg.rpcTimeResolution > 0.);
0259               const double stripTime =
0260                   (*Digitization::Gauss{calibCfg.rpcTimeResolution}(hit.time(),
0261                                                                     rndEngine))
0262                       .first;
0263               newSp.setTime(stripTime);
0264             }
0265             newSp.setCovariance(
0266                 calibCfg.rpcPhiStripPitch, calibCfg.rpcEtaStripPitch,
0267                 m_cfg.digitizeTime ? calibCfg.rpcTimeResolution : 0.);
0268
0269             break;
0270           }
0271           /// Endcap strips not yet available
0272           case SurfaceBounds::BoundsType::eTrapezoid:
0273             break;
0274           default:
0275             convertSp = false;
0276         }
0277         /// Implement a dead time
0278         if (convertSp) {
0279           newSp.defineCoordinates(
0280               Vector3{parentTrf * smearedHit},
0281               Vector3{parentTrf.linear() * Vector3::UnitY()},
0282               Vector3{parentTrf.linear() * Vector3::UnitX()});
0283           MuonId_t id{};
0284           /// @todo Refine me using the volume name
0285           id.setChamber(MuonId_t::StationName::BIS,
0286                         hit.position().z() > 0 ? MuonId_t::DetSide::A
0287                                                : MuonId_t::DetSide::C,
0288                         1, MuonId_t::TechField::Rpc);
0289           id.setCoordFlags(true, true);
0290           newSp.setId(id);
0291         }
0292
0293         break;
0294       }
0295       case Straw: {
0296         auto closeApproach =
0297             lineIntersect<3>(Vector3::Zero(), Vector3::UnitZ(), locPos, locDir);
0298         const auto nominalPos = closeApproach.position();
0299         const double unsmearedR = fastHypot(nominalPos.x(), nominalPos.y());
0300         ACTS_VERBOSE("Hit is recorded in a straw detector, R: "
0301                      << unsmearedR << ", " << bounds);
0302
0303         const double uncert = calibrator().driftRadiusUncert(unsmearedR);
0304         /// Reject unsmearable hits
0305         if (uncert <= std::numeric_limits<double>::epsilon()) {
0306           convertSp = false;
0307           break;
0308         }
0309         const double driftR =
0310             (*Digitization::Gauss{uncert}(unsmearedR, rndEngine)).first;
0311         // bounds
0312         const auto& lBounds = static_cast<const LineBounds&>(bounds);
0313         const double maxR = lBounds.get(LineBounds::eR);
0314         const double maxZ = lBounds.get(LineBounds::eHalfLengthZ);
0315         /// The generated hit is unphysical
0316         if (driftR < 0. || driftR > maxR || std::abs(nominalPos.z()) > maxZ) {
0317           convertSp = false;
0318           break;
0319         }
0320         if (auto insertItr =
0321                 strawTimes.insert(std::make_pair(hitId, hit.time()));
0322             !insertItr.second) {
0323           if (hit.time() - insertItr.first->second > m_cfg.strawDeadTime) {
0324             insertItr.first->second = hit.time();
0325           } else {
0326             convertSp = false;
0327             break;
0328           }
0329         }
0330
0331         newSp.setRadius(driftR);
0332         newSp.setCovariance(square(0.5 * maxZ),
0333                             calibrator().driftRadiusUncert(driftR), 0.);
0334         newSp.defineCoordinates(Vector3{parentTrf.translation()},
0335                                 Vector3{parentTrf.linear() * Vector3::UnitZ()},
0336                                 Vector3{parentTrf.linear() * Vector3::UnitX()});
0337         MuonId_t id{};
0338         /// @todo Refine me using the volume name
0339         id.setChamber(MuonId_t::StationName::BIS,
0340                       hit.position().z() > 0 ? MuonId_t::DetSide::A
0341                                              : MuonId_t::DetSide::C,
0342                       1, MuonId_t::TechField::Mdt);
0343         id.setCoordFlags(true, false);
0344         newSp.setId(id);
0345         break;
0346       }
0347       ///
0348       default:
0349         convertSp = false;
0350     }
0351
0352     if (convertSp) {
0353       spacePointsPerChamber[toChamberId(hitId)].push_back(std::move(newSp));
0354     }
0355   }
0356
0357   for (auto& [volId, bucket] : spacePointsPerChamber) {
0358     std::ranges::sort(bucket,
0359                       [](const MuonSpacePoint& a, const MuonSpacePoint& b) {
0360                         return a.localPosition().z() < b.localPosition().z();
0361                       });
0362     if (logger().doPrint(Logging::Level::VERBOSE)) {
0363       std::stringstream sstr{};
0364       for (const auto& spacePoint : bucket) {
0365         sstr << " *** " << spacePoint << std::endl;
0366       }
0367       ACTS_VERBOSE("Safe " << bucket.size() << " space points for chamber "
0368                            << volId << "\n"
0369                            << sstr.str());
0370     }
0371     visualizeBucket(ctx, gctx, bucket);
0372     outSpacePoints.push_back(std::move(bucket));
0373   }
0374   m_outputSpacePoints(ctx, std::move(outSpacePoints));
0375
0376   return ProcessCode::SUCCESS;
0377 }
0378
0379 RangeXD<2, double> MuonSpacePointDigitizer::canvasRanges(
0380     const MuonSpacePointBucket& bucket) const {
0381   RangeXD<2, double> ranges{
0382       filledArray<double, 2>(std::numeric_limits<double>::max()),
0383       filledArray<double, 2>(-std::numeric_limits<double>::max())};
0384   // Extra margin such that the canvas axes don't overlap with the depicted
0385   // measurements
0386   constexpr double extra = 3._cm;
0387   for (const auto& sp : bucket) {
0388     const Vector3& pos = sp.localPosition();
0389     ranges.expand(0, pos.z() - extra, pos.z() + extra);
0390     ranges.expand(1, pos.y() - extra, pos.y() + extra);
0391   }
0392   return ranges;
0393 }
0394
0395 bool MuonSpacePointDigitizer::isSurfaceToDraw(
0396     const Acts::GeometryContext& gctx, const Surface& surface,
0397     const RangeXD<2, double>& canvasBoundaries) const {
0398   // Draw only active surfaces
0399   if (surface.associatedDetectorElement() == nullptr) {
0400     return false;
0401   }
0402   // surface position in the frame
0403   const Vector3 pos =
0404       toSpacePointFrame(gctx, surface.geometryId()).translation();
0405   const auto& bounds = surface.bounds();
0406
0407   if (surface.type() == Surface::SurfaceType::Plane) {
0408     const double hL{halfHeight(bounds)};
0409     // check whether the surface is inside the visible range
0410     const double minZ = std::max(pos.z() - hL, canvasBoundaries.min(0));
0411     const double maxZ = std::min(pos.z() + hL, canvasBoundaries.max(0));
0412     // The maximum is below the left side of the strip plane
0413     // or the minimum is above the right side of the plane
0414     return maxZ > pos.z() - hL && minZ < pos.z() + hL &&
0415            pos.y() > canvasBoundaries.min(1) &&
0416            pos.y() < canvasBoundaries.max(1);
0417   } else if (surface.type() == Surface::SurfaceType::Straw) {
0418     const double r = static_cast<const LineBounds&>(bounds).get(LineBounds::eR);
0419     // Check that the straw surface is well embedded on the canvas
0420     return pos.y() - r > canvasBoundaries.min(1) &&
0421            pos.y() + r < canvasBoundaries.max(1) &&
0422            pos.z() - r > canvasBoundaries.min(0) &&
0423            pos.z() + r < canvasBoundaries.max(0);
0424   }
0425
0426   return false;
0427 }
0428 void MuonSpacePointDigitizer::visualizeBucket(
0429     const AlgorithmContext& ctx, const GeometryContext& gctx,
0430     const MuonSpacePointBucket& bucket) const {
0431   if (!m_cfg.dumpVisualization) {
0432     return;
0433   }
0434   auto canvas = std::make_unique<TCanvas>("can", "can", 600, 600);
0435   canvas->cd();
0436
0437   const GeometryIdentifier chambId = toChamberId(bucket.front().geometryId());
0438
0439   std::vector<std::unique_ptr<TObject>> primitives{};
0440
0441   const RangeXD<2, double> canvasBound{canvasRanges(bucket)};
0442   /// Draw the frame
0443   auto frame = std::make_unique<TH2I>("frame", "frame;z [mm];y [mm]", 1,
0444                                       canvasBound.min(0), canvasBound.max(0), 1,
0445                                       canvasBound.min(1), canvasBound.max(1));
0446   frame->Draw("AXIS");
0447
0448   // Loop over all surfaces inside the chamber volume to draw the ones covered
0449   // by the canvas
0450   const TrackingVolume* chambVolume = trackingGeometry().findVolume(chambId);
0451   assert(chambVolume != nullptr);
0452   chambVolume->apply(overloaded{
0453       [this, &canvasBound, &gctx, &primitives](const Surface& surface) {
0454         if (!isSurfaceToDraw(gctx, surface, canvasBound)) {
0455           return;
0456         }
0457         const Vector3 pos =
0458             toSpacePointFrame(gctx, surface.geometryId()).translation();
0459         const auto& bounds = surface.bounds();
0460         if (surface.type() == Surface::SurfaceType::Plane) {
0461           const double hL{halfHeight(bounds)};
0462           const double minZ = std::max(pos.z() - hL, canvasBound.min(0));
0463           const double maxZ = std::min(pos.z() + hL, canvasBound.max(0));
0464           primitives.push_back(drawBox(0.5 * (minZ + maxZ), pos.y(),
0465                                        maxZ - minZ, 0.3_cm, kBlack, 0));
0466
0467         } else if (surface.type() == Surface::SurfaceType::Straw) {
0468           const double r =
0469               static_cast<const LineBounds&>(bounds).get(LineBounds::eR);
0470           primitives.push_back(drawCircle(pos.z(), pos.y(), r, kBlack, 0));
0471         }
0472       }});
0473
0474   for (auto& sp : bucket) {
0475     const Vector3& pos = sp.localPosition();
0476     if (sp.isStraw()) {
0477       primitives.push_back(
0478           drawCircle(pos.z(), pos.y(), sp.driftRadius(), kRed, 0));
0479     } else {
0480       primitives.push_back(
0481           drawBox(pos.z(), pos.y(), 3._cm, 0.5_cm, kRed, 1001));
0482     }
0483   }
0484   // Finally draw the muon trajectory
0485   const SimHitContainer& gotSimHits = m_inputSimHits(ctx);
0486   const SimParticleContainer& simParticles = m_inputParticles(ctx);
0487   for (const auto& simHit : gotSimHits) {
0488     if (chambId != toChamberId(simHit.geometryId())) {
0489       continue;
0490     }
0491     const auto simPartItr = simParticles.find(simHit.particleId());
0492     if (simPartItr == simParticles.end() ||
0493         (*simPartItr).hypothesis() != ParticleHypothesis::muon()) {
0494       continue;
0495     }
0496     const auto toSpTrf = toSpacePointFrame(gctx, simHit.geometryId()) *
0497                          trackingGeometry()
0498                              .findSurface(simHit.geometryId())
0499                              ->transform(gctx)
0500                              .inverse();
0501     const Vector3 pos = toSpTrf * simHit.position();
0502     const Vector3 dir = toSpTrf.linear() * simHit.direction();
0503     constexpr double arrowL = 1._cm;
0504     const Vector3 start = pos - 0.5 * arrowL * dir;
0505     const Vector3 end = pos + 0.5 * arrowL * dir;
0506     auto arrow =
0507         std::make_unique<TArrow>(start.z(), start.y(), end.z(), end.y(), 0.03);
0508     arrow->SetLineColor(kBlue + 1);
0509     arrow->SetLineWidth(1);
0510     arrow->SetLineStyle(kSolid);
0511     primitives.push_back(std::move(arrow));
0512   }
0513
0514   for (auto& prim : primitives) {
0515     prim->Draw();
0516   }
0517   canvas->SaveAs(
0518       std::format("Event_{}_{}.pdf", ctx.eventNumber, chambVolume->volumeName())
0519           .c_str());
0520 }
0521
0522 }  // namespace ActsExamples