src/EventData/MuonSpacePointCalibrator.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/EventData/MuonSpacePointCalibrator.hpp"
0010
0011 #include "Acts/Surfaces/detail/LineHelper.hpp"
0012 #include "Acts/Utilities/detail/Polynomials.hpp"
0013
0014 namespace ActsExamples {
0015
0016 using TechField = MuonSpacePoint::MuonId::TechField;
0017
0018 MuonSpacePointCalibrator::MuonSpacePointCalibrator(
0019     const Config& cfg, std::unique_ptr<const Acts::Logger> logger)
0020     : m_logger{std::move(logger)}, m_cfg{cfg} {}
0021 MuonSpacePointCalibrator::CalibSpCont_t MuonSpacePointCalibrator::calibrate(
0022     const Acts::CalibrationContext& ctx, const Acts::Vector3& trackPos,
0023     const Acts::Vector3& trackDir, const double trackT0,
0024     const UnCalibSpVec_t& uncalibCont) const {
0025   CalibSpCont_t outContainer{};
0026   outContainer.reserve(uncalibCont.size());
0027   for (const MuonSpacePoint* calibMe : uncalibCont) {
0028     calibrate(ctx, trackPos, trackDir, trackT0, *calibMe, outContainer);
0029   }
0030   return outContainer;
0031 }
0032
0033 void MuonSpacePointCalibrator::calibrate(
0034     const Acts::CalibrationContext& /*ctx*/, const Acts::Vector3& trackPos,
0035     const Acts::Vector3& trackDir, const double trackT0,
0036     const MuonSpacePoint& spacePoint, CalibSpCont_t& outContainer) const {
0037   using namespace Acts::detail::LineHelper;
0038   auto calibSp = std::make_unique<MuonSpacePoint>(spacePoint);
0039   const Acts::Intersection3D closePoint =
0040       lineIntersect<3>(trackPos, trackDir, spacePoint.localPosition(),
0041                        spacePoint.sensorDirection());
0042
0043   switch (spacePoint.id().technology()) {
0044     using enum TechField;
0045     case Mdt: {
0046       if (closePoint.isValid()) {
0047         calibSp->defineCoordinates(Acts::Vector3{closePoint.position()},
0048                                    Acts::Vector3{spacePoint.sensorDirection()},
0049                                    Acts::Vector3{spacePoint.toNextSensor()});
0050
0051         double driftT = driftTime(spacePoint.driftRadius());
0052         if (m_cfg.includeTrackT0) {
0053           driftT -= trackT0;
0054         }
0055         calibSp->setRadius(driftRadius(driftT));
0056       }
0057       break;
0058     }
0059     /** Micromega & sTGC technology don't re-adjust the time */
0060     case Tgc:
0061     case Mm:
0062     case sTgc: {
0063       if (closePoint.isValid() &&
0064           (!spacePoint.id().measuresEta() || !spacePoint.id().measuresPhi())) {
0065         calibSp->defineCoordinates(Acts::Vector3{closePoint.position()},
0066                                    Acts::Vector3{spacePoint.sensorDirection()},
0067                                    Acts::Vector3{spacePoint.toNextSensor()});
0068       }
0069       break;
0070     }
0071     case Rpc: {
0072       /** TODO: Add time adjustment taking the signal propagation time into
0073        * account */
0074       if (closePoint.isValid() &&
0075           (!spacePoint.id().measuresEta() || !spacePoint.id().measuresPhi())) {
0076         calibSp->defineCoordinates(Acts::Vector3{closePoint.position()},
0077                                    Acts::Vector3{spacePoint.sensorDirection()},
0078                                    Acts::Vector3{spacePoint.toNextSensor()});
0079       }
0080       break;
0081     }
0082     default:
0083       break;
0084   }
0085   outContainer.push_back(std::move(calibSp));
0086 }
0087
0088 std::optional<double> MuonSpacePointCalibrator::expandPolySeries(
0089     const double x, unsigned derivative, const CalibPolyType polyType,
0090     const double upperBound, const double lowerBound,
0091     const std::vector<double>& coeffs) const {
0092   if (x < lowerBound || x > upperBound) {
0093     ACTS_VERBOSE("Parsed x:" << x << " is outside of the configured range ["
0094                              << lowerBound << ", " << upperBound << "]");
0095     return std::nullopt;
0096   }
0097   const double intervalNorm = 2. / (upperBound - lowerBound);
0098   const double xNorm = intervalNorm * (x - 0.5 * (upperBound + lowerBound));
0099   const double chainFactor = Acts::pow(intervalNorm, derivative);
0100   double result{0.};
0101   const auto N = static_cast<std::uint32_t>(coeffs.size());
0102   for (std::uint32_t k = 0; k < N; ++k) {
0103     switch (polyType) {
0104       case CalibPolyType::Chebychev:
0105         result +=
0106             Acts::detail::chebychevPolyTn(xNorm, k, derivative) * coeffs[k];
0107         break;
0108       case CalibPolyType::Legendre:
0109         result += Acts::detail::legendrePoly(xNorm, k, derivative) * coeffs[k];
0110         break;
0111     }
0112   }
0113   return result * chainFactor;
0114 }
0115
0116 double MuonSpacePointCalibrator::driftRadius(const double driftTime) const {
0117   return expandPolySeries(driftTime, 0u, m_cfg.rtPolyType, m_cfg.maxDriftT,
0118                           m_cfg.minDriftT, m_cfg.rtCoefficients)
0119       .value_or(0.);
0120 }
0121 double MuonSpacePointCalibrator::driftRadiusUncert(
0122     const double driftRadius) const {
0123   return expandPolySeries(driftTime(driftRadius), 0u, m_cfg.rtPolyType,
0124                           m_cfg.maxDriftT, m_cfg.minDriftT,
0125                           m_cfg.rtCoefficients)
0126       .value_or(0.);
0127 }
0128 double MuonSpacePointCalibrator::driftVelocity(const double driftTime) const {
0129   return expandPolySeries(driftTime, 1u, m_cfg.rtPolyType, m_cfg.maxDriftT,
0130                           m_cfg.minDriftT, m_cfg.rtCoefficients)
0131       .value_or(0.);
0132 }
0133 double MuonSpacePointCalibrator::driftAcceleration(
0134     const double driftTime) const {
0135   return expandPolySeries(driftTime, 2u, m_cfg.rtPolyType, m_cfg.maxDriftT,
0136                           m_cfg.minDriftT, m_cfg.rtCoefficients)
0137       .value_or(m_cfg.minTubeR);
0138 }
0139 double MuonSpacePointCalibrator::driftTime(const double driftRadius) const {
0140   return expandPolySeries(driftRadius, 0, m_cfg.trPolyType, m_cfg.maxTubeR,
0141                           m_cfg.minTubeR, m_cfg.trCoefficients)
0142       .value_or(m_cfg.minDriftT);
0143 }
0144
0145 double MuonSpacePointCalibrator::driftVelocity(
0146     const Acts::CalibrationContext& /*ctx*/, const MuonSpacePoint& sp) const {
0147   return sp.id().technology() == MuonSpacePoint::MuonId::TechField::Mdt
0148              ? driftVelocity(driftTime(sp.driftRadius()))
0149              : 0.;
0150 }
0151
0152 double MuonSpacePointCalibrator::driftAcceleration(
0153     const Acts::CalibrationContext& /*ctx*/, const MuonSpacePoint& sp) const {
0154   return sp.id().technology() == MuonSpacePoint::MuonId::TechField::Mdt
0155              ? driftAcceleration(driftTime(sp.driftRadius()))
0156              : 0.;
0157 }
0158
0159 }  // namespace ActsExamples