src/Validation/ResPlotTool.cpp

0001 // This file is part of the Acts project.
0002 //
0003 // Copyright (C) 2019-2023 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 http://mozilla.org/MPL/2.0/.
0008
0009 #include "ActsExamples/Validation/ResPlotTool.hpp"
0010
0011 #include "Acts/Definitions/TrackParametrization.hpp"
0012 #include "Acts/Surfaces/Surface.hpp"
0013 #include "Acts/Utilities/Helpers.hpp"
0014 #include "Acts/Utilities/Result.hpp"
0015 #include "ActsFatras/EventData/Particle.hpp"
0016
0017 #include <algorithm>
0018 #include <cmath>
0019 #include <optional>
0020 #include <ostream>
0021
0022 #include <TH1.h>
0023 #include <TH2.h>
0024 #include <TString.h>
0025
0026 ActsExamples::ResPlotTool::ResPlotTool(
0027     const ActsExamples::ResPlotTool::Config& cfg, Acts::Logging::Level lvl)
0028     : m_cfg(cfg), m_logger(Acts::getDefaultLogger("ResPlotTool", lvl)) {}
0029
0030 void ActsExamples::ResPlotTool::book(
0031     ResPlotTool::ResPlotCache& resPlotCache) const {
0032   PlotHelpers::Binning bEta = m_cfg.varBinning.at("Eta");
0033   PlotHelpers::Binning bPt = m_cfg.varBinning.at("Pt");
0034   PlotHelpers::Binning bPull = m_cfg.varBinning.at("Pull");
0035
0036   ACTS_DEBUG("Initialize the histograms for residual and pull plots");
0037   for (unsigned int parID = 0; parID < Acts::eBoundSize; parID++) {
0038     std::string parName = m_cfg.paramNames.at(parID);
0039
0040     std::string parResidual = "Residual_" + parName;
0041     // Binning for residual is parameter dependent
0042     PlotHelpers::Binning bResidual = m_cfg.varBinning.at(parResidual);
0043
0044     // residual distributions
0045     resPlotCache.res[parName] = PlotHelpers::bookHisto(
0046         Form("res_%s", parName.c_str()),
0047         Form("Residual of %s", parName.c_str()), bResidual);
0048     // residual vs eta scatter plots
0049     resPlotCache.res_vs_eta[parName] = PlotHelpers::bookHisto(
0050         Form("res_%s_vs_eta", parName.c_str()),
0051         Form("Residual of %s vs eta", parName.c_str()), bEta, bResidual);
0052     // residual mean in each eta bin
0053     resPlotCache.resMean_vs_eta[parName] = PlotHelpers::bookHisto(
0054         Form("resmean_%s_vs_eta", parName.c_str()),
0055         Form("Residual mean of %s", parName.c_str()), bEta);
0056     // residual width in each eta bin
0057     resPlotCache.resWidth_vs_eta[parName] = PlotHelpers::bookHisto(
0058         Form("reswidth_%s_vs_eta", parName.c_str()),
0059         Form("Residual width of %s", parName.c_str()), bEta);
0060     // residual vs pT scatter plots
0061     resPlotCache.res_vs_pT[parName] = PlotHelpers::bookHisto(
0062         Form("res_%s_vs_pT", parName.c_str()),
0063         Form("Residual of %s vs pT", parName.c_str()), bPt, bResidual);
0064     // residual mean in each pT bin
0065     resPlotCache.resMean_vs_pT[parName] = PlotHelpers::bookHisto(
0066         Form("resmean_%s_vs_pT", parName.c_str()),
0067         Form("Residual mean of %s", parName.c_str()), bPt);
0068     // residual width in each pT bin
0069     resPlotCache.resWidth_vs_pT[parName] = PlotHelpers::bookHisto(
0070         Form("reswidth_%s_vs_pT", parName.c_str()),
0071         Form("Residual width of %s", parName.c_str()), bPt);
0072
0073     // pull distritutions
0074     resPlotCache.pull[parName] =
0075         PlotHelpers::bookHisto(Form("pull_%s", parName.c_str()),
0076                                Form("Pull of %s", parName.c_str()), bPull);
0077     // pull vs eta scatter plots
0078     resPlotCache.pull_vs_eta[parName] = PlotHelpers::bookHisto(
0079         Form("pull_%s_vs_eta", parName.c_str()),
0080         Form("Pull of %s vs eta", parName.c_str()), bEta, bPull);
0081     // pull mean in each eta bin
0082     resPlotCache.pullMean_vs_eta[parName] =
0083         PlotHelpers::bookHisto(Form("pullmean_%s_vs_eta", parName.c_str()),
0084                                Form("Pull mean of %s", parName.c_str()), bEta);
0085     // pull width in each eta bin
0086     resPlotCache.pullWidth_vs_eta[parName] =
0087         PlotHelpers::bookHisto(Form("pullwidth_%s_vs_eta", parName.c_str()),
0088                                Form("Pull width of %s", parName.c_str()), bEta);
0089     // pull vs pT scatter plots
0090     resPlotCache.pull_vs_pT[parName] = PlotHelpers::bookHisto(
0091         Form("pull_%s_vs_pT", parName.c_str()),
0092         Form("Pull of %s vs pT", parName.c_str()), bPt, bPull);
0093     // pull mean in each pT bin
0094     resPlotCache.pullMean_vs_pT[parName] =
0095         PlotHelpers::bookHisto(Form("pullmean_%s_vs_pT", parName.c_str()),
0096                                Form("Pull mean of %s", parName.c_str()), bPt);
0097     // pull width in each pT bin
0098     resPlotCache.pullWidth_vs_pT[parName] =
0099         PlotHelpers::bookHisto(Form("pullwidth_%s_vs_pT", parName.c_str()),
0100                                Form("Pull width of %s", parName.c_str()), bPt);
0101   }
0102 }
0103
0104 void ActsExamples::ResPlotTool::clear(ResPlotCache& resPlotCache) const {
0105   ACTS_DEBUG("Delete the hists.");
0106   for (unsigned int parID = 0; parID < Acts::eBoundSize; parID++) {
0107     std::string parName = m_cfg.paramNames.at(parID);
0108     delete resPlotCache.res.at(parName);
0109     delete resPlotCache.res_vs_eta.at(parName);
0110     delete resPlotCache.resMean_vs_eta.at(parName);
0111     delete resPlotCache.resWidth_vs_eta.at(parName);
0112     delete resPlotCache.res_vs_pT.at(parName);
0113     delete resPlotCache.resMean_vs_pT.at(parName);
0114     delete resPlotCache.resWidth_vs_pT.at(parName);
0115     delete resPlotCache.pull.at(parName);
0116     delete resPlotCache.pull_vs_eta.at(parName);
0117     delete resPlotCache.pullMean_vs_eta.at(parName);
0118     delete resPlotCache.pullWidth_vs_eta.at(parName);
0119     delete resPlotCache.pull_vs_pT.at(parName);
0120     delete resPlotCache.pullMean_vs_pT.at(parName);
0121     delete resPlotCache.pullWidth_vs_pT.at(parName);
0122   }
0123 }
0124
0125 void ActsExamples::ResPlotTool::write(
0126     const ResPlotTool::ResPlotCache& resPlotCache) const {
0127   ACTS_DEBUG("Write the hists to output file.");
0128   for (unsigned int parID = 0; parID < Acts::eBoundSize; parID++) {
0129     std::string parName = m_cfg.paramNames.at(parID);
0130     resPlotCache.res.at(parName)->Write();
0131     resPlotCache.res_vs_eta.at(parName)->Write();
0132     resPlotCache.resMean_vs_eta.at(parName)->Write();
0133     resPlotCache.resWidth_vs_eta.at(parName)->Write();
0134     resPlotCache.res_vs_pT.at(parName)->Write();
0135     resPlotCache.resMean_vs_pT.at(parName)->Write();
0136     resPlotCache.resWidth_vs_pT.at(parName)->Write();
0137     resPlotCache.pull.at(parName)->Write();
0138     resPlotCache.pull_vs_eta.at(parName)->Write();
0139     resPlotCache.pullMean_vs_eta.at(parName)->Write();
0140     resPlotCache.pullWidth_vs_eta.at(parName)->Write();
0141     resPlotCache.pull_vs_pT.at(parName)->Write();
0142     resPlotCache.pullMean_vs_pT.at(parName)->Write();
0143     resPlotCache.pullWidth_vs_pT.at(parName)->Write();
0144   }
0145 }
0146
0147 void ActsExamples::ResPlotTool::fill(
0148     ResPlotTool::ResPlotCache& resPlotCache, const Acts::GeometryContext& gctx,
0149     const ActsFatras::Particle& truthParticle,
0150     const Acts::BoundTrackParameters& fittedParamters) const {
0151   using ParametersVector = Acts::BoundTrackParameters::ParametersVector;
0152   using Acts::VectorHelpers::eta;
0153   using Acts::VectorHelpers::perp;
0154   using Acts::VectorHelpers::phi;
0155   using Acts::VectorHelpers::theta;
0156
0157   // get the fitted parameter (at perigee surface) and its error
0158   auto trackParameter = fittedParamters.parameters();
0159
0160   // get the perigee surface
0161   auto pSurface = &fittedParamters.referenceSurface();
0162
0163   // get the truth position and momentum
0164   ParametersVector truthParameter = ParametersVector::Zero();
0165
0166   // get the truth perigee parameter
0167   auto lpResult = pSurface->globalToLocal(gctx, truthParticle.position(),
0168                                           truthParticle.direction());
0169   if (lpResult.ok()) {
0170     truthParameter[Acts::BoundIndices::eBoundLoc0] =
0171         lpResult.value()[Acts::BoundIndices::eBoundLoc0];
0172     truthParameter[Acts::BoundIndices::eBoundLoc1] =
0173         lpResult.value()[Acts::BoundIndices::eBoundLoc1];
0174   } else {
0175     ACTS_ERROR("Global to local transformation did not succeed.");
0176   }
0177   truthParameter[Acts::BoundIndices::eBoundPhi] =
0178       phi(truthParticle.direction());
0179   truthParameter[Acts::BoundIndices::eBoundTheta] =
0180       theta(truthParticle.direction());
0181   truthParameter[Acts::BoundIndices::eBoundQOverP] =
0182       truthParticle.charge() / truthParticle.absoluteMomentum();
0183   truthParameter[Acts::BoundIndices::eBoundTime] = truthParticle.time();
0184
0185   // get the truth eta and pT
0186   const auto truthEta = eta(truthParticle.direction());
0187   const auto truthPt = truthParticle.transverseMomentum();
0188
0189   // fill the histograms for residual and pull
0190   for (unsigned int parID = 0; parID < Acts::eBoundSize; parID++) {
0191     std::string parName = m_cfg.paramNames.at(parID);
0192     float residual = trackParameter[parID] - truthParameter[parID];
0193     PlotHelpers::fillHisto(resPlotCache.res.at(parName), residual);
0194     PlotHelpers::fillHisto(resPlotCache.res_vs_eta.at(parName), truthEta,
0195                            residual);
0196     PlotHelpers::fillHisto(resPlotCache.res_vs_pT.at(parName), truthPt,
0197                            residual);
0198
0199     if (fittedParamters.covariance().has_value()) {
0200       auto covariance = *fittedParamters.covariance();
0201       if (covariance(parID, parID) > 0) {
0202         float pull = residual / sqrt(covariance(parID, parID));
0203         PlotHelpers::fillHisto(resPlotCache.pull[parName], pull);
0204         PlotHelpers::fillHisto(resPlotCache.pull_vs_eta.at(parName), truthEta,
0205                                pull);
0206         PlotHelpers::fillHisto(resPlotCache.pull_vs_pT.at(parName), truthPt,
0207                                pull);
0208       } else {
0209         ACTS_WARNING("Fitted track parameter :" << parName
0210                                                 << " has negative covariance = "
0211                                                 << covariance(parID, parID));
0212       }
0213     } else {
0214       ACTS_WARNING("Fitted track parameter :" << parName
0215                                               << " has no covariance");
0216     }
0217   }
0218 }
0219
0220 // get the mean and width of residual/pull in each eta/pT bin and fill them into
0221 // histograms
0222 void ActsExamples::ResPlotTool::refinement(
0223     ResPlotTool::ResPlotCache& resPlotCache) const {
0224   PlotHelpers::Binning bEta = m_cfg.varBinning.at("Eta");
0225   PlotHelpers::Binning bPt = m_cfg.varBinning.at("Pt");
0226   for (unsigned int parID = 0; parID < Acts::eBoundSize; parID++) {
0227     std::string parName = m_cfg.paramNames.at(parID);
0228     // refine the plots vs eta
0229     for (int j = 1; j <= static_cast<int>(bEta.nBins()); j++) {
0230       TH1D* temp_res = resPlotCache.res_vs_eta.at(parName)->ProjectionY(
0231           Form("%s_projy_bin%d", "Residual_vs_eta_Histo", j), j, j);
0232       PlotHelpers::anaHisto(temp_res, j,
0233                             resPlotCache.resMean_vs_eta.at(parName),
0234                             resPlotCache.resWidth_vs_eta.at(parName));
0235
0236       TH1D* temp_pull = resPlotCache.pull_vs_eta.at(parName)->ProjectionY(
0237           Form("%s_projy_bin%d", "Pull_vs_eta_Histo", j), j, j);
0238       PlotHelpers::anaHisto(temp_pull, j,
0239                             resPlotCache.pullMean_vs_eta.at(parName),
0240                             resPlotCache.pullWidth_vs_eta.at(parName));
0241     }
0242
0243     // refine the plots vs pT
0244     for (int j = 1; j <= static_cast<int>(bPt.nBins()); j++) {
0245       TH1D* temp_res = resPlotCache.res_vs_pT.at(parName)->ProjectionY(
0246           Form("%s_projy_bin%d", "Residual_vs_pT_Histo", j), j, j);
0247       PlotHelpers::anaHisto(temp_res, j, resPlotCache.resMean_vs_pT.at(parName),
0248                             resPlotCache.resWidth_vs_pT.at(parName));
0249
0250       TH1D* temp_pull = resPlotCache.pull_vs_pT.at(parName)->ProjectionY(
0251           Form("%s_projy_bin%d", "Pull_vs_pT_Histo", j), j, j);
0252       PlotHelpers::anaHisto(temp_pull, j,
0253                             resPlotCache.pullMean_vs_pT.at(parName),
0254                             resPlotCache.pullWidth_vs_pT.at(parName));
0255     }
0256   }
0257 }