src/Validation/ResPlotTool.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/Validation/ResPlotTool.hpp"
0010
0011 #include "Acts/Definitions/Algebra.hpp"
0012 #include "Acts/Definitions/TrackParametrization.hpp"
0013 #include "Acts/Surfaces/Surface.hpp"
0014 #include "Acts/Utilities/Intersection.hpp"
0015 #include "Acts/Utilities/Result.hpp"
0016 #include "Acts/Utilities/detail/periodic.hpp"
0017
0018 #include <cmath>
0019 #include <cstdint>
0020 #include <format>
0021 #include <numbers>
0022 #include <stdexcept>
0023
0024 namespace ActsExamples {
0025
0026 static constexpr double nan = std::numeric_limits<double>::quiet_NaN();
0027
0028 /// Reduce a residual of one bound parameter to its minimal equivalent value,
0029 /// which only matters for the periodic phi.
0030 static double wrapResidual(std::size_t index, double residual) {
0031   if (index != Acts::eBoundPhi) {
0032     return residual;
0033   }
0034   return Acts::detail::difference_periodic(residual, 0., 2 * std::numbers::pi);
0035 }
0036
0037 ResPlotTool::ResPlotTool(const ResPlotTool::Config& cfg,
0038                          Acts::Logging::Level lvl)
0039     : m_cfg(cfg), m_logger(Acts::getDefaultLogger("ResPlotTool", lvl)) {
0040   // `varBinning.at` would only report the key type, not the missing key
0041   const auto binning = [this](const std::string& key) -> const AxisVariant& {
0042     const auto it = m_cfg.varBinning.find(key);
0043     if (it == m_cfg.varBinning.end()) {
0044       throw std::invalid_argument("ResPlotTool: missing binning for '" + key +
0045                                   "'");
0046     }
0047     return it->second;
0048   };
0049
0050   const auto& etaAxis = binning("Eta");
0051   const auto& phiAxis = binning("Phi");
0052   const auto& ptAxis = binning("Pt");
0053   const auto& pullAxis = binning("Pull");
0054
0055   if (m_cfg.paramNames.size() != Acts::eBoundSize) {
0056     throw std::invalid_argument(
0057         "ResPlotTool: expected one name per bound parameter, the caller has to "
0058         "fill in the empty default");
0059   }
0060
0061   ACTS_DEBUG("Initialize the histograms for residual and pull plots");
0062
0063   std::vector<std::string> allParamNames = m_cfg.paramNames;
0064   allParamNames.push_back(m_cfg.qOverPtName);
0065   allParamNames.push_back(m_cfg.relQoverPtName);
0066
0067   for (const std::string& parName : allParamNames) {
0068     const auto& residualAxis = binning("Residual_" + parName);
0069
0070     // residual distributions
0071     m_res.emplace(parName, Acts::Experimental::Histogram1(
0072                                std::format("res_{}", parName),
0073                                std::format("Residual of {}", parName),
0074                                std::array{residualAxis}));
0075
0076     // residual vs eta scatter plots
0077     m_resVsEta.emplace(parName,
0078                        Acts::Experimental::Histogram2(
0079                            std::format("res_{}_vs_eta", parName),
0080                            std::format("Residual of {} vs eta", parName),
0081                            std::array{etaAxis, residualAxis}));
0082
0083     // residual vs pT scatter plots
0084     m_resVsPt.emplace(parName, Acts::Experimental::Histogram2(
0085                                    std::format("res_{}_vs_pT", parName),
0086                                    std::format("Residual of {} vs pT", parName),
0087                                    std::array{ptAxis, residualAxis}));
0088
0089     // residual vs eta-phi scatter plots
0090     m_resVsEtaPhi.emplace(parName,
0091                           Acts::Experimental::Histogram3(
0092                               std::format("res_{}_vs_eta_phi", parName),
0093                               std::format("Residual of {} vs eta-phi", parName),
0094                               std::array{etaAxis, phiAxis, residualAxis}));
0095
0096     // residual vs eta-pT scatter plots
0097     m_resVsEtaPt.emplace(parName,
0098                          Acts::Experimental::Histogram3(
0099                              std::format("res_{}_vs_eta_pT", parName),
0100                              std::format("Residual of {} vs eta-pT", parName),
0101                              std::array{etaAxis, ptAxis, residualAxis}));
0102
0103     // pull distributions
0104     m_pull.emplace(
0105         parName, Acts::Experimental::Histogram1(
0106                      std::format("pull_{}", parName),
0107                      std::format("Pull of {}", parName), std::array{pullAxis}));
0108
0109     // pull vs eta scatter plots
0110     m_pullVsEta.emplace(parName, Acts::Experimental::Histogram2(
0111                                      std::format("pull_{}_vs_eta", parName),
0112                                      std::format("Pull of {} vs eta", parName),
0113                                      std::array{etaAxis, pullAxis}));
0114
0115     // pull vs pT scatter plots
0116     m_pullVsPt.emplace(parName, Acts::Experimental::Histogram2(
0117                                     std::format("pull_{}_vs_pT", parName),
0118                                     std::format("Pull of {} vs pT", parName),
0119                                     std::array{ptAxis, pullAxis}));
0120
0121     // pull vs eta-phi scatter plots
0122     m_pullVsEtaPhi.emplace(parName,
0123                            Acts::Experimental::Histogram3(
0124                                std::format("pull_{}_vs_eta_phi", parName),
0125                                std::format("Pull of {} vs eta-phi", parName),
0126                                std::array{etaAxis, phiAxis, pullAxis}));
0127
0128     // pull vs eta-pT scatter plots
0129     m_pullVsEtaPt.emplace(parName,
0130                           Acts::Experimental::Histogram3(
0131                               std::format("pull_{}_vs_eta_pT", parName),
0132                               std::format("Pull of {} vs eta-pT", parName),
0133                               std::array{etaAxis, ptAxis, pullAxis}));
0134   }
0135 }
0136
0137 void ResPlotTool::fill(const Acts::GeometryContext& gctx,
0138                        const SimParticleState& truthParticle,
0139                        const Acts::BoundTrackParameters& fittedParamters) {
0140   using Acts::VectorHelpers::eta;
0141   using Acts::VectorHelpers::perp;
0142   using Acts::VectorHelpers::phi;
0143   using Acts::VectorHelpers::theta;
0144
0145   using enum Acts::BoundIndices;
0146
0147   // get the perigee surface
0148   const Acts::Surface& pSurface = fittedParamters.referenceSurface();
0149
0150   // get the truth parameter at the perigee surface
0151   Acts::BoundVector truthParameters = Acts::BoundVector::Zero();
0152   const Acts::Intersection3D intersection =
0153       pSurface
0154           .intersect(gctx, truthParticle.position(), truthParticle.direction())
0155           .closest();
0156   if (intersection.isValid()) {
0157     const Acts::Result<Acts::Vector2> lpResult = pSurface.globalToLocal(
0158         gctx, intersection.position(), truthParticle.direction());
0159     assert(lpResult.ok());
0160
0161     truthParameters[eBoundLoc0] = lpResult.value()[eBoundLoc0];
0162     truthParameters[eBoundLoc1] = lpResult.value()[eBoundLoc1];
0163   } else {
0164     ACTS_ERROR("Cannot get the truth perigee parameter");
0165   }
0166   truthParameters[eBoundPhi] = phi(truthParticle.direction());
0167   truthParameters[eBoundTheta] = theta(truthParticle.direction());
0168   truthParameters[eBoundQOverP] = truthParticle.qOverP();
0169   truthParameters[eBoundTime] = truthParticle.time();
0170
0171   // bin on the particle, not the bound parameters, which would round-trip the
0172   // direction through phi/theta
0173   fill(truthParameters,
0174        Binning{eta(truthParticle.direction()), phi(truthParticle.direction()),
0175                truthParticle.transverseMomentum()},
0176        truthParticle.charge(), truthParticle.absoluteCharge(), fittedParamters);
0177 }
0178
0179 void ResPlotTool::fill(const Acts::BoundTrackParameters& truthParameters,
0180                        const Acts::BoundTrackParameters& fittedParameters) {
0181   using Acts::VectorHelpers::eta;
0182   using Acts::VectorHelpers::phi;
0183
0184   if (truthParameters.referenceSurface() !=
0185       fittedParameters.referenceSurface()) {
0186     throw std::invalid_argument(
0187         "ResPlotTool: truth and fitted parameters are expressed on different "
0188         "reference surfaces");
0189   }
0190
0191   const double truthCharge = truthParameters.charge();
0192   fill(truthParameters.parameters(),
0193        Binning{eta(truthParameters.direction()),
0194                phi(truthParameters.direction()),
0195                truthParameters.transverseMomentum()},
0196        truthCharge, std::abs(truthCharge), fittedParameters);
0197 }
0198
0199 void ResPlotTool::fill(const Acts::BoundVector& truthVector,
0200                        const Binning& binning, double truthCharge,
0201                        double truthAbsCharge,
0202                        const Acts::BoundTrackParameters& fittedParameters) {
0203   using enum Acts::BoundIndices;
0204
0205   const double truthEta = binning.eta;
0206   const double truthPhi = binning.phi;
0207   const double truthPt = binning.pt;
0208
0209   // get the fitted parameter and its error
0210   const Acts::BoundVector& trackParameters = fittedParameters.parameters();
0211   const Acts::BoundMatrix& trackCovariance =
0212       fittedParameters.covariance().value_or(Acts::BoundMatrix::Zero());
0213
0214   // fill the histograms for residual and pull
0215   for (unsigned int paramId = 0; paramId < Acts::eBoundSize; paramId++) {
0216     const std::string& parName = m_cfg.paramNames.at(paramId);
0217
0218     const double residual =
0219         wrapResidual(paramId, trackParameters[paramId] - truthVector[paramId]);
0220     fillResidual(parName, residual, truthEta, truthPhi, truthPt);
0221
0222     const double var = trackCovariance(paramId, paramId);
0223
0224     const double pull = var > 0 ? residual / std::sqrt(var) : nan;
0225     fillPull(parName, pull, truthEta, truthPhi, truthPt);
0226   }
0227
0228   // `reco(q/pT)` and `true(pT/q) * reco(q/pT)` residual and pull
0229   {
0230     const double truthQoverPt = truthCharge / truthPt;
0231     const double truthPtOverAbsQ = truthPt / truthAbsCharge;
0232     const double recoQoverPt =
0233         trackParameters[eBoundQOverP] / std::sin(trackParameters[eBoundTheta]);
0234     const double residualQoverPt = recoQoverPt - truthQoverPt;
0235     fillResidual(m_cfg.qOverPtName, residualQoverPt, truthEta, truthPhi,
0236                  truthPt);
0237
0238     const double residualRelQoverPt = truthPtOverAbsQ * residualQoverPt;
0239     fillResidual(m_cfg.relQoverPtName, residualRelQoverPt, truthEta, truthPhi,
0240                  truthPt);
0241
0242     const double covarianceQoverPt = [&]() -> double {
0243       const Acts::Vector2 jacobian{
0244           -recoQoverPt / std::tan(trackParameters[eBoundTheta]),
0245           1 / std::sin(trackParameters[eBoundTheta])};
0246       const Acts::SquareMatrix2 covariance = trackCovariance(
0247           {eBoundTheta, eBoundQOverP}, {eBoundTheta, eBoundQOverP});
0248       return jacobian.transpose() * covariance * jacobian;
0249     }();
0250     const double covarianceRelQoverPt =
0251         Acts::square(truthPtOverAbsQ) * covarianceQoverPt;
0252
0253     const double pullQoverPt =
0254         covarianceQoverPt > 0 ? residualQoverPt / std::sqrt(covarianceQoverPt)
0255                               : nan;
0256     fillPull(m_cfg.qOverPtName, pullQoverPt, truthEta, truthPhi, truthPt);
0257
0258     const double pullRelQoverPt =
0259         covarianceRelQoverPt > 0
0260             ? residualRelQoverPt / std::sqrt(covarianceRelQoverPt)
0261             : nan;
0262     fillPull(m_cfg.relQoverPtName, pullRelQoverPt, truthEta, truthPhi, truthPt);
0263   }
0264 }
0265
0266 void ResPlotTool::fill(const Binning& binning,
0267                        const Acts::VariableBoundSubspaceHelper& subspace,
0268                        const Acts::BoundVector& residuals,
0269                        const Acts::BoundMatrix& residualCovariance) {
0270   for (const std::uint8_t index : subspace) {
0271     const std::string& parName = m_cfg.paramNames.at(index);
0272
0273     const double residual = wrapResidual(index, residuals[index]);
0274     fillResidual(parName, residual, binning.eta, binning.phi, binning.pt);
0275
0276     // `V - HPH^T` is not positive definite, so there is not always a pull
0277     const double var = residualCovariance(index, index);
0278     const double pull = var > 0 ? residual / std::sqrt(var) : nan;
0279     fillPull(parName, pull, binning.eta, binning.phi, binning.pt);
0280   }
0281 }
0282
0283 void ResPlotTool::fillResidual(const std::string& paramName, double residual,
0284                                double truthEta, double truthPhi,
0285                                double truthPt) {
0286   m_res.at(paramName).fill({residual});
0287   m_resVsEta.at(paramName).fill({truthEta, residual});
0288   m_resVsPt.at(paramName).fill({truthPt, residual});
0289   m_resVsEtaPhi.at(paramName).fill({truthEta, truthPhi, residual});
0290   m_resVsEtaPt.at(paramName).fill({truthEta, truthPt, residual});
0291 }
0292
0293 void ResPlotTool::fillPull(const std::string& paramName, double pull,
0294                            double truthEta, double truthPhi, double truthPt) {
0295   m_pull.at(paramName).fill({pull});
0296   m_pullVsEta.at(paramName).fill({truthEta, pull});
0297   m_pullVsPt.at(paramName).fill({truthPt, pull});
0298   m_pullVsEtaPhi.at(paramName).fill({truthEta, truthPhi, pull});
0299   m_pullVsEtaPt.at(paramName).fill({truthEta, truthPt, pull});
0300 }
0301
0302 }  // namespace ActsExamples