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