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
0017 #include <cmath>
0018 #include <format>
0019
0020 namespace ActsExamples {
0021
0022 static constexpr double nan = std::numeric_limits<double>::quiet_NaN();
0023
0024 ResPlotTool::ResPlotTool(const ResPlotTool::Config& cfg,
0025                          Acts::Logging::Level lvl)
0026     : m_cfg(cfg), m_logger(Acts::getDefaultLogger("ResPlotTool", lvl)) {
0027   const auto& etaAxis = m_cfg.varBinning.at("Eta");
0028   const auto& phiAxis = m_cfg.varBinning.at("Phi");
0029   const auto& ptAxis = m_cfg.varBinning.at("Pt");
0030   const auto& pullAxis = m_cfg.varBinning.at("Pull");
0031
0032   ACTS_DEBUG("Initialize the histograms for residual and pull plots");
0033
0034   std::vector<std::string> allParamNames = m_cfg.paramNames;
0035   allParamNames.push_back(m_cfg.qOverPtName);
0036   allParamNames.push_back(m_cfg.relQoverPtName);
0037
0038   for (const std::string& parName : allParamNames) {
0039     const std::string parResidual = "Residual_" + parName;
0040     const auto& residualAxis = m_cfg.varBinning.at(parResidual);
0041
0042     // residual distributions
0043     m_res.emplace(parName, Acts::Experimental::Histogram1(
0044                                std::format("res_{}", parName),
0045                                std::format("Residual of {}", parName),
0046                                std::array{residualAxis}));
0047
0048     // residual vs eta scatter plots
0049     m_resVsEta.emplace(parName,
0050                        Acts::Experimental::Histogram2(
0051                            std::format("res_{}_vs_eta", parName),
0052                            std::format("Residual of {} vs eta", parName),
0053                            std::array{etaAxis, residualAxis}));
0054
0055     // residual vs pT scatter plots
0056     m_resVsPt.emplace(parName, Acts::Experimental::Histogram2(
0057                                    std::format("res_{}_vs_pT", parName),
0058                                    std::format("Residual of {} vs pT", parName),
0059                                    std::array{ptAxis, residualAxis}));
0060
0061     // residual vs eta-phi scatter plots
0062     m_resVsEtaPhi.emplace(parName,
0063                           Acts::Experimental::Histogram3(
0064                               std::format("res_{}_vs_eta_phi", parName),
0065                               std::format("Residual of {} vs eta-phi", parName),
0066                               std::array{etaAxis, phiAxis, residualAxis}));
0067
0068     // residual vs eta-pT scatter plots
0069     m_resVsEtaPt.emplace(parName,
0070                          Acts::Experimental::Histogram3(
0071                              std::format("res_{}_vs_eta_pT", parName),
0072                              std::format("Residual of {} vs eta-pT", parName),
0073                              std::array{etaAxis, ptAxis, residualAxis}));
0074
0075     // pull distributions
0076     m_pull.emplace(
0077         parName, Acts::Experimental::Histogram1(
0078                      std::format("pull_{}", parName),
0079                      std::format("Pull of {}", parName), std::array{pullAxis}));
0080
0081     // pull vs eta scatter plots
0082     m_pullVsEta.emplace(parName, Acts::Experimental::Histogram2(
0083                                      std::format("pull_{}_vs_eta", parName),
0084                                      std::format("Pull of {} vs eta", parName),
0085                                      std::array{etaAxis, pullAxis}));
0086
0087     // pull vs pT scatter plots
0088     m_pullVsPt.emplace(parName, Acts::Experimental::Histogram2(
0089                                     std::format("pull_{}_vs_pT", parName),
0090                                     std::format("Pull of {} vs pT", parName),
0091                                     std::array{ptAxis, pullAxis}));
0092
0093     // pull vs eta-phi scatter plots
0094     m_pullVsEtaPhi.emplace(parName,
0095                            Acts::Experimental::Histogram3(
0096                                std::format("pull_{}_vs_eta_phi", parName),
0097                                std::format("Pull of {} vs eta-phi", parName),
0098                                std::array{etaAxis, phiAxis, pullAxis}));
0099
0100     // pull vs eta-pT scatter plots
0101     m_pullVsEtaPt.emplace(parName,
0102                           Acts::Experimental::Histogram3(
0103                               std::format("pull_{}_vs_eta_pT", parName),
0104                               std::format("Pull of {} vs eta-pT", parName),
0105                               std::array{etaAxis, ptAxis, pullAxis}));
0106   }
0107 }
0108
0109 void ResPlotTool::fill(const Acts::GeometryContext& gctx,
0110                        const SimParticleState& truthParticle,
0111                        const Acts::BoundTrackParameters& fittedParamters) {
0112   using Acts::VectorHelpers::eta;
0113   using Acts::VectorHelpers::perp;
0114   using Acts::VectorHelpers::phi;
0115   using Acts::VectorHelpers::theta;
0116
0117   using enum Acts::BoundIndices;
0118
0119   // get the fitted parameter (at perigee surface) and its error
0120   const Acts::BoundVector& trackParameters = fittedParamters.parameters();
0121   const Acts::BoundMatrix& trackCovariance =
0122       fittedParamters.covariance().value_or(Acts::BoundMatrix::Zero());
0123
0124   // get the perigee surface
0125   const Acts::Surface& pSurface = fittedParamters.referenceSurface();
0126
0127   // get the truth parameter at the perigee surface
0128   Acts::BoundVector truthParameters = Acts::BoundVector::Zero();
0129   const Acts::Intersection3D intersection =
0130       pSurface
0131           .intersect(gctx, truthParticle.position(), truthParticle.direction())
0132           .closest();
0133   if (intersection.isValid()) {
0134     const Acts::Result<Acts::Vector2> lpResult = pSurface.globalToLocal(
0135         gctx, intersection.position(), truthParticle.direction());
0136     assert(lpResult.ok());
0137
0138     truthParameters[eBoundLoc0] = lpResult.value()[eBoundLoc0];
0139     truthParameters[eBoundLoc1] = lpResult.value()[eBoundLoc1];
0140   } else {
0141     ACTS_ERROR("Cannot get the truth perigee parameter");
0142   }
0143   truthParameters[eBoundPhi] = phi(truthParticle.direction());
0144   truthParameters[eBoundTheta] = theta(truthParticle.direction());
0145   truthParameters[eBoundQOverP] = truthParticle.qOverP();
0146   truthParameters[eBoundTime] = truthParticle.time();
0147
0148   // get the truth eta and pT
0149   const double truthEta = eta(truthParticle.direction());
0150   const double truthPhi = phi(truthParticle.direction());
0151   const double truthPt = truthParticle.transverseMomentum();
0152
0153   // fill the histograms for residual and pull
0154   for (unsigned int paramId = 0; paramId < Acts::eBoundSize; paramId++) {
0155     const std::string& parName = m_cfg.paramNames.at(paramId);
0156
0157     const double residual = trackParameters[paramId] - truthParameters[paramId];
0158     fillResidual(parName, residual, truthEta, truthPhi, truthPt);
0159
0160     const double var = trackCovariance(paramId, paramId);
0161
0162     const double pull = var > 0 ? residual / std::sqrt(var) : nan;
0163     fillPull(parName, pull, truthEta, truthPhi, truthPt);
0164   }
0165
0166   // `reco(q/pT)` and `true(pT/q) * reco(q/pT)` residual and pull
0167   {
0168     const double truthQoverPt = truthParticle.charge() / truthPt;
0169     const double truthPtOverQ = truthPt / truthParticle.charge();
0170     const double recoQoverPt =
0171         trackParameters[eBoundQOverP] / std::sin(trackParameters[eBoundTheta]);
0172     const double residualQoverPt = recoQoverPt - truthQoverPt;
0173     fillResidual(m_cfg.qOverPtName, residualQoverPt, truthEta, truthPhi,
0174                  truthPt);
0175
0176     const double residualRelQoverPt = truthPtOverQ * residualQoverPt;
0177     fillResidual(m_cfg.relQoverPtName, residualRelQoverPt, truthEta, truthPhi,
0178                  truthPt);
0179
0180     const double covarianceQoverPt = [&]() -> double {
0181       const Acts::Vector2 jacobian{
0182           -recoQoverPt / std::tan(trackParameters[eBoundTheta]),
0183           1 / std::sin(trackParameters[eBoundTheta])};
0184       const Acts::SquareMatrix2 covariance = trackCovariance(
0185           {eBoundTheta, eBoundQOverP}, {eBoundTheta, eBoundQOverP});
0186       return jacobian.transpose() * covariance * jacobian;
0187     }();
0188     const double covarianceRelQoverPt =
0189         Acts::square(truthPtOverQ) * covarianceQoverPt;
0190
0191     const double pullQoverPt =
0192         covarianceQoverPt > 0 ? residualQoverPt / std::sqrt(covarianceQoverPt)
0193                               : nan;
0194     fillPull(m_cfg.qOverPtName, pullQoverPt, truthEta, truthPhi, truthPt);
0195
0196     const double pullRelQoverPt =
0197         covarianceRelQoverPt > 0
0198             ? residualRelQoverPt / std::sqrt(covarianceRelQoverPt)
0199             : nan;
0200     fillPull(m_cfg.relQoverPtName, pullRelQoverPt, truthEta, truthPhi, truthPt);
0201   }
0202 }
0203
0204 void ResPlotTool::fillResidual(const std::string& paramName, double residual,
0205                                double truthEta, double truthPhi,
0206                                double truthPt) {
0207   m_res.at(paramName).fill({residual});
0208   m_resVsEta.at(paramName).fill({truthEta, residual});
0209   m_resVsPt.at(paramName).fill({truthPt, residual});
0210   m_resVsEtaPhi.at(paramName).fill({truthEta, truthPhi, residual});
0211   m_resVsEtaPt.at(paramName).fill({truthEta, truthPt, residual});
0212 }
0213
0214 void ResPlotTool::fillPull(const std::string& paramName, double pull,
0215                            double truthEta, double truthPhi, double truthPt) {
0216   m_pull.at(paramName).fill({pull});
0217   m_pullVsEta.at(paramName).fill({truthEta, pull});
0218   m_pullVsPt.at(paramName).fill({truthPt, pull});
0219   m_pullVsEtaPhi.at(paramName).fill({truthEta, truthPhi, pull});
0220   m_pullVsEtaPt.at(paramName).fill({truthEta, truthPt, pull});
0221 }
0222
0223 }  // namespace ActsExamples