Back to home page

EIC code displayed by LXR
 
 

    

File indexing completed on 2025-10-26 07:55:29 

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/Intersection.hpp"
0014 #include "Acts/Utilities/Result.hpp"
0015 #include "ActsExamples/EventData/SimParticle.hpp"
0016 
0017 #include <cmath>
0018 #include <optional>
0019 #include <ostream>
0020 
0021 #include <TH1.h>
0022 #include <TH2.h>
0023 #include <TString.h>
0024 
0025 namespace ActsExamples {
0026 
0027 ResPlotTool::ResPlotTool(const ResPlotTool::Config& cfg,
0028                          Acts::Logging::Level lvl)
0029     : m_cfg(cfg), m_logger(Acts::getDefaultLogger("ResPlotTool", lvl)) {}
0030 
0031 void ResPlotTool::book(Cache& cache) 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     cache.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     cache.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     cache.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     cache.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     cache.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     cache.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     cache.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     cache.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     cache.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     cache.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     cache.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     cache.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     cache.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     cache.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 ResPlotTool::clear(Cache& cache) 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 cache.res.at(parName);
0109     delete cache.res_vs_eta.at(parName);
0110     delete cache.resMean_vs_eta.at(parName);
0111     delete cache.resWidth_vs_eta.at(parName);
0112     delete cache.res_vs_pT.at(parName);
0113     delete cache.resMean_vs_pT.at(parName);
0114     delete cache.resWidth_vs_pT.at(parName);
0115     delete cache.pull.at(parName);
0116     delete cache.pull_vs_eta.at(parName);
0117     delete cache.pullMean_vs_eta.at(parName);
0118     delete cache.pullWidth_vs_eta.at(parName);
0119     delete cache.pull_vs_pT.at(parName);
0120     delete cache.pullMean_vs_pT.at(parName);
0121     delete cache.pullWidth_vs_pT.at(parName);
0122   }
0123 }
0124 
0125 void ResPlotTool::write(const Cache& cache) const {
0126   ACTS_DEBUG("Write the hists to output file.");
0127   for (unsigned int parID = 0; parID < Acts::eBoundSize; parID++) {
0128     std::string parName = m_cfg.paramNames.at(parID);
0129     cache.res.at(parName)->Write();
0130     cache.res_vs_eta.at(parName)->Write();
0131     cache.resMean_vs_eta.at(parName)->Write();
0132     cache.resWidth_vs_eta.at(parName)->Write();
0133     cache.res_vs_pT.at(parName)->Write();
0134     cache.resMean_vs_pT.at(parName)->Write();
0135     cache.resWidth_vs_pT.at(parName)->Write();
0136     cache.pull.at(parName)->Write();
0137     cache.pull_vs_eta.at(parName)->Write();
0138     cache.pullMean_vs_eta.at(parName)->Write();
0139     cache.pullWidth_vs_eta.at(parName)->Write();
0140     cache.pull_vs_pT.at(parName)->Write();
0141     cache.pullMean_vs_pT.at(parName)->Write();
0142     cache.pullWidth_vs_pT.at(parName)->Write();
0143   }
0144 }
0145 
0146 void ResPlotTool::fill(
0147     Cache& cache, const Acts::GeometryContext& gctx,
0148     const SimParticleState& truthParticle,
0149     const Acts::BoundTrackParameters& fittedParamters) const {
0150   using ParametersVector = Acts::BoundTrackParameters::ParametersVector;
0151   using Acts::VectorHelpers::eta;
0152   using Acts::VectorHelpers::perp;
0153   using Acts::VectorHelpers::phi;
0154   using Acts::VectorHelpers::theta;
0155 
0156   // get the fitted parameter (at perigee surface) and its error
0157   Acts::BoundVector trackParameter = fittedParamters.parameters();
0158 
0159   // get the perigee surface
0160   const Acts::Surface& pSurface = fittedParamters.referenceSurface();
0161 
0162   // get the truth position and momentum
0163   ParametersVector truthParameter = ParametersVector::Zero();
0164 
0165   // get the truth perigee parameter
0166   Acts::Intersection3D intersection =
0167       pSurface
0168           .intersect(gctx, truthParticle.position(), truthParticle.direction())
0169           .closest();
0170   if (intersection.isValid()) {
0171     auto lpResult = pSurface.globalToLocal(gctx, intersection.position(),
0172                                            truthParticle.direction());
0173     assert(lpResult.ok());
0174 
0175     truthParameter[Acts::BoundIndices::eBoundLoc0] =
0176         lpResult.value()[Acts::BoundIndices::eBoundLoc0];
0177     truthParameter[Acts::BoundIndices::eBoundLoc1] =
0178         lpResult.value()[Acts::BoundIndices::eBoundLoc1];
0179   } else {
0180     ACTS_ERROR("Cannot get the truth perigee parameter");
0181   }
0182   truthParameter[Acts::BoundIndices::eBoundPhi] =
0183       phi(truthParticle.direction());
0184   truthParameter[Acts::BoundIndices::eBoundTheta] =
0185       theta(truthParticle.direction());
0186   truthParameter[Acts::BoundIndices::eBoundQOverP] = truthParticle.qOverP();
0187   truthParameter[Acts::BoundIndices::eBoundTime] = truthParticle.time();
0188 
0189   // get the truth eta and pT
0190   const auto truthEta = eta(truthParticle.direction());
0191   const auto truthPt = truthParticle.transverseMomentum();
0192 
0193   // fill the histograms for residual and pull
0194   for (unsigned int parID = 0; parID < Acts::eBoundSize; parID++) {
0195     std::string parName = m_cfg.paramNames.at(parID);
0196     float residual = trackParameter[parID] - truthParameter[parID];
0197     PlotHelpers::fillHisto(cache.res.at(parName), residual);
0198     PlotHelpers::fillHisto(cache.res_vs_eta.at(parName), truthEta, residual);
0199     PlotHelpers::fillHisto(cache.res_vs_pT.at(parName), truthPt, residual);
0200 
0201     if (fittedParamters.covariance().has_value()) {
0202       auto covariance = *fittedParamters.covariance();
0203       if (covariance(parID, parID) > 0) {
0204         float pull = residual / sqrt(covariance(parID, parID));
0205         PlotHelpers::fillHisto(cache.pull[parName], pull);
0206         PlotHelpers::fillHisto(cache.pull_vs_eta.at(parName), truthEta, pull);
0207         PlotHelpers::fillHisto(cache.pull_vs_pT.at(parName), truthPt, 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 ResPlotTool::refinement(Cache& cache) const {
0223   PlotHelpers::Binning bEta = m_cfg.varBinning.at("Eta");
0224   PlotHelpers::Binning bPt = m_cfg.varBinning.at("Pt");
0225   for (unsigned int parID = 0; parID < Acts::eBoundSize; parID++) {
0226     std::string parName = m_cfg.paramNames.at(parID);
0227     // refine the plots vs eta
0228     for (int j = 1; j <= static_cast<int>(bEta.nBins()); j++) {
0229       TH1D* temp_res = cache.res_vs_eta.at(parName)->ProjectionY(
0230           Form("%s_projy_bin%d", "Residual_vs_eta_Histo", j), j, j);
0231       PlotHelpers::anaHisto(temp_res, j, cache.resMean_vs_eta.at(parName),
0232                             cache.resWidth_vs_eta.at(parName));
0233 
0234       TH1D* temp_pull = cache.pull_vs_eta.at(parName)->ProjectionY(
0235           Form("%s_projy_bin%d", "Pull_vs_eta_Histo", j), j, j);
0236       PlotHelpers::anaHisto(temp_pull, j, cache.pullMean_vs_eta.at(parName),
0237                             cache.pullWidth_vs_eta.at(parName));
0238     }
0239 
0240     // refine the plots vs pT
0241     for (int j = 1; j <= static_cast<int>(bPt.nBins()); j++) {
0242       TH1D* temp_res = cache.res_vs_pT.at(parName)->ProjectionY(
0243           Form("%s_projy_bin%d", "Residual_vs_pT_Histo", j), j, j);
0244       PlotHelpers::anaHisto(temp_res, j, cache.resMean_vs_pT.at(parName),
0245                             cache.resWidth_vs_pT.at(parName));
0246 
0247       TH1D* temp_pull = cache.pull_vs_pT.at(parName)->ProjectionY(
0248           Form("%s_projy_bin%d", "Pull_vs_pT_Histo", j), j, j);
0249       PlotHelpers::anaHisto(temp_pull, j, cache.pullMean_vs_pT.at(parName),
0250                             cache.pullWidth_vs_pT.at(parName));
0251     }
0252   }
0253 }
0254 
0255 }  // namespace ActsExamples
[ Source navigation ] [ Diff markup ] [ Identifier search ] [ general search ]

This page was automatically generated by the 2.3.7 LXR engine.

The LXR team
Valid CSS 2.1! Valid HTML 4.01!