Back to home page

EIC code displayed by LXR
 
 

    

File indexing completed on 2025-01-18 09:11:49 

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 ActsExamples::ResPlotTool::ResPlotTool(
0025     const ActsExamples::ResPlotTool::Config& cfg, Acts::Logging::Level lvl)
0026     : m_cfg(cfg), m_logger(Acts::getDefaultLogger("ResPlotTool", lvl)) {}
0027 
0028 void ActsExamples::ResPlotTool::book(
0029     ResPlotTool::ResPlotCache& resPlotCache) const {
0030   PlotHelpers::Binning bEta = m_cfg.varBinning.at("Eta");
0031   PlotHelpers::Binning bPt = m_cfg.varBinning.at("Pt");
0032   PlotHelpers::Binning bPull = m_cfg.varBinning.at("Pull");
0033 
0034   ACTS_DEBUG("Initialize the histograms for residual and pull plots");
0035   for (unsigned int parID = 0; parID < Acts::eBoundSize; parID++) {
0036     std::string parName = m_cfg.paramNames.at(parID);
0037 
0038     std::string parResidual = "Residual_" + parName;
0039     // Binning for residual is parameter dependent
0040     PlotHelpers::Binning bResidual = m_cfg.varBinning.at(parResidual);
0041 
0042     // residual distributions
0043     resPlotCache.res[parName] = PlotHelpers::bookHisto(
0044         Form("res_%s", parName.c_str()),
0045         Form("Residual of %s", parName.c_str()), bResidual);
0046     // residual vs eta scatter plots
0047     resPlotCache.res_vs_eta[parName] = PlotHelpers::bookHisto(
0048         Form("res_%s_vs_eta", parName.c_str()),
0049         Form("Residual of %s vs eta", parName.c_str()), bEta, bResidual);
0050     // residual mean in each eta bin
0051     resPlotCache.resMean_vs_eta[parName] = PlotHelpers::bookHisto(
0052         Form("resmean_%s_vs_eta", parName.c_str()),
0053         Form("Residual mean of %s", parName.c_str()), bEta);
0054     // residual width in each eta bin
0055     resPlotCache.resWidth_vs_eta[parName] = PlotHelpers::bookHisto(
0056         Form("reswidth_%s_vs_eta", parName.c_str()),
0057         Form("Residual width of %s", parName.c_str()), bEta);
0058     // residual vs pT scatter plots
0059     resPlotCache.res_vs_pT[parName] = PlotHelpers::bookHisto(
0060         Form("res_%s_vs_pT", parName.c_str()),
0061         Form("Residual of %s vs pT", parName.c_str()), bPt, bResidual);
0062     // residual mean in each pT bin
0063     resPlotCache.resMean_vs_pT[parName] = PlotHelpers::bookHisto(
0064         Form("resmean_%s_vs_pT", parName.c_str()),
0065         Form("Residual mean of %s", parName.c_str()), bPt);
0066     // residual width in each pT bin
0067     resPlotCache.resWidth_vs_pT[parName] = PlotHelpers::bookHisto(
0068         Form("reswidth_%s_vs_pT", parName.c_str()),
0069         Form("Residual width of %s", parName.c_str()), bPt);
0070 
0071     // pull distritutions
0072     resPlotCache.pull[parName] =
0073         PlotHelpers::bookHisto(Form("pull_%s", parName.c_str()),
0074                                Form("Pull of %s", parName.c_str()), bPull);
0075     // pull vs eta scatter plots
0076     resPlotCache.pull_vs_eta[parName] = PlotHelpers::bookHisto(
0077         Form("pull_%s_vs_eta", parName.c_str()),
0078         Form("Pull of %s vs eta", parName.c_str()), bEta, bPull);
0079     // pull mean in each eta bin
0080     resPlotCache.pullMean_vs_eta[parName] =
0081         PlotHelpers::bookHisto(Form("pullmean_%s_vs_eta", parName.c_str()),
0082                                Form("Pull mean of %s", parName.c_str()), bEta);
0083     // pull width in each eta bin
0084     resPlotCache.pullWidth_vs_eta[parName] =
0085         PlotHelpers::bookHisto(Form("pullwidth_%s_vs_eta", parName.c_str()),
0086                                Form("Pull width of %s", parName.c_str()), bEta);
0087     // pull vs pT scatter plots
0088     resPlotCache.pull_vs_pT[parName] = PlotHelpers::bookHisto(
0089         Form("pull_%s_vs_pT", parName.c_str()),
0090         Form("Pull of %s vs pT", parName.c_str()), bPt, bPull);
0091     // pull mean in each pT bin
0092     resPlotCache.pullMean_vs_pT[parName] =
0093         PlotHelpers::bookHisto(Form("pullmean_%s_vs_pT", parName.c_str()),
0094                                Form("Pull mean of %s", parName.c_str()), bPt);
0095     // pull width in each pT bin
0096     resPlotCache.pullWidth_vs_pT[parName] =
0097         PlotHelpers::bookHisto(Form("pullwidth_%s_vs_pT", parName.c_str()),
0098                                Form("Pull width of %s", parName.c_str()), bPt);
0099   }
0100 }
0101 
0102 void ActsExamples::ResPlotTool::clear(ResPlotCache& resPlotCache) const {
0103   ACTS_DEBUG("Delete the hists.");
0104   for (unsigned int parID = 0; parID < Acts::eBoundSize; parID++) {
0105     std::string parName = m_cfg.paramNames.at(parID);
0106     delete resPlotCache.res.at(parName);
0107     delete resPlotCache.res_vs_eta.at(parName);
0108     delete resPlotCache.resMean_vs_eta.at(parName);
0109     delete resPlotCache.resWidth_vs_eta.at(parName);
0110     delete resPlotCache.res_vs_pT.at(parName);
0111     delete resPlotCache.resMean_vs_pT.at(parName);
0112     delete resPlotCache.resWidth_vs_pT.at(parName);
0113     delete resPlotCache.pull.at(parName);
0114     delete resPlotCache.pull_vs_eta.at(parName);
0115     delete resPlotCache.pullMean_vs_eta.at(parName);
0116     delete resPlotCache.pullWidth_vs_eta.at(parName);
0117     delete resPlotCache.pull_vs_pT.at(parName);
0118     delete resPlotCache.pullMean_vs_pT.at(parName);
0119     delete resPlotCache.pullWidth_vs_pT.at(parName);
0120   }
0121 }
0122 
0123 void ActsExamples::ResPlotTool::write(
0124     const ResPlotTool::ResPlotCache& resPlotCache) 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     resPlotCache.res.at(parName)->Write();
0129     resPlotCache.res_vs_eta.at(parName)->Write();
0130     resPlotCache.resMean_vs_eta.at(parName)->Write();
0131     resPlotCache.resWidth_vs_eta.at(parName)->Write();
0132     resPlotCache.res_vs_pT.at(parName)->Write();
0133     resPlotCache.resMean_vs_pT.at(parName)->Write();
0134     resPlotCache.resWidth_vs_pT.at(parName)->Write();
0135     resPlotCache.pull.at(parName)->Write();
0136     resPlotCache.pull_vs_eta.at(parName)->Write();
0137     resPlotCache.pullMean_vs_eta.at(parName)->Write();
0138     resPlotCache.pullWidth_vs_eta.at(parName)->Write();
0139     resPlotCache.pull_vs_pT.at(parName)->Write();
0140     resPlotCache.pullMean_vs_pT.at(parName)->Write();
0141     resPlotCache.pullWidth_vs_pT.at(parName)->Write();
0142   }
0143 }
0144 
0145 void ActsExamples::ResPlotTool::fill(
0146     ResPlotTool::ResPlotCache& resPlotCache, 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(resPlotCache.res.at(parName), residual);
0197     PlotHelpers::fillHisto(resPlotCache.res_vs_eta.at(parName), truthEta,
0198                            residual);
0199     PlotHelpers::fillHisto(resPlotCache.res_vs_pT.at(parName), truthPt,
0200                            residual);
0201 
0202     if (fittedParamters.covariance().has_value()) {
0203       auto covariance = *fittedParamters.covariance();
0204       if (covariance(parID, parID) > 0) {
0205         float pull = residual / sqrt(covariance(parID, parID));
0206         PlotHelpers::fillHisto(resPlotCache.pull[parName], pull);
0207         PlotHelpers::fillHisto(resPlotCache.pull_vs_eta.at(parName), truthEta,
0208                                pull);
0209         PlotHelpers::fillHisto(resPlotCache.pull_vs_pT.at(parName), truthPt,
0210                                pull);
0211       } else {
0212         ACTS_WARNING("Fitted track parameter :" << parName
0213                                                 << " has negative covariance = "
0214                                                 << covariance(parID, parID));
0215       }
0216     } else {
0217       ACTS_WARNING("Fitted track parameter :" << parName
0218                                               << " has no covariance");
0219     }
0220   }
0221 }
0222 
0223 // get the mean and width of residual/pull in each eta/pT bin and fill them into
0224 // histograms
0225 void ActsExamples::ResPlotTool::refinement(
0226     ResPlotTool::ResPlotCache& resPlotCache) const {
0227   PlotHelpers::Binning bEta = m_cfg.varBinning.at("Eta");
0228   PlotHelpers::Binning bPt = m_cfg.varBinning.at("Pt");
0229   for (unsigned int parID = 0; parID < Acts::eBoundSize; parID++) {
0230     std::string parName = m_cfg.paramNames.at(parID);
0231     // refine the plots vs eta
0232     for (int j = 1; j <= static_cast<int>(bEta.nBins()); j++) {
0233       TH1D* temp_res = resPlotCache.res_vs_eta.at(parName)->ProjectionY(
0234           Form("%s_projy_bin%d", "Residual_vs_eta_Histo", j), j, j);
0235       PlotHelpers::anaHisto(temp_res, j,
0236                             resPlotCache.resMean_vs_eta.at(parName),
0237                             resPlotCache.resWidth_vs_eta.at(parName));
0238 
0239       TH1D* temp_pull = resPlotCache.pull_vs_eta.at(parName)->ProjectionY(
0240           Form("%s_projy_bin%d", "Pull_vs_eta_Histo", j), j, j);
0241       PlotHelpers::anaHisto(temp_pull, j,
0242                             resPlotCache.pullMean_vs_eta.at(parName),
0243                             resPlotCache.pullWidth_vs_eta.at(parName));
0244     }
0245 
0246     // refine the plots vs pT
0247     for (int j = 1; j <= static_cast<int>(bPt.nBins()); j++) {
0248       TH1D* temp_res = resPlotCache.res_vs_pT.at(parName)->ProjectionY(
0249           Form("%s_projy_bin%d", "Residual_vs_pT_Histo", j), j, j);
0250       PlotHelpers::anaHisto(temp_res, j, resPlotCache.resMean_vs_pT.at(parName),
0251                             resPlotCache.resWidth_vs_pT.at(parName));
0252 
0253       TH1D* temp_pull = resPlotCache.pull_vs_pT.at(parName)->ProjectionY(
0254           Form("%s_projy_bin%d", "Pull_vs_pT_Histo", j), j, j);
0255       PlotHelpers::anaHisto(temp_pull, j,
0256                             resPlotCache.pullMean_vs_pT.at(parName),
0257                             resPlotCache.pullWidth_vs_pT.at(parName));
0258     }
0259   }
0260 }
[ 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!