File indexing completed on 2026-07-16 08:08:01
0001
0002
0003
0004
0005
0006
0007
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
0047
0048
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
0056
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
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
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
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 }
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
0153
0154 const TrackingVolume* volume =
0155 trackingGeometry().findVolume(toChamberId(hitId));
0156 assert(volume != nullptr);
0157
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
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
0187
0188
0189
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
0199 const Transform3 parentTrf{toSpacePointFrame(gctx, moduleGeoId)};
0200
0201 const auto& bounds = hitSurf->bounds();
0202
0203
0204 for (auto h = moduleSimHits.begin(); h != moduleSimHits.end(); ++h) {
0205 const auto& simHit = *h;
0206
0207
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
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
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
0265
0266
0267 stripTimes.insert(std::make_pair(
0268 moduleGeoId, std::array{smearedHit[ePos0], smearedHit[ePos1],
0269 simHit.time()}));
0270
0271
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
0287 case SurfaceBounds::BoundsType::eTrapezoid:
0288 break;
0289 default:
0290 convertSp = false;
0291 }
0292
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
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
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
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
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
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
0409
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
0423 if (!surface.isSensitive()) {
0424 return false;
0425 }
0426
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
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
0437
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
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
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
0473
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
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 }