src/Validation/EffPlotTool.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/EffPlotTool.hpp"
0010
0011 #include "Acts/Utilities/VectorHelpers.hpp"
0012 #include "Acts/Utilities/Zip.hpp"
0013 #include "ActsExamples/EventData/SimParticle.hpp"
0014
0015 #include <format>
0016 #include <limits>
0017
0018 using Acts::VectorHelpers::eta;
0019 using Acts::VectorHelpers::perp;
0020 using Acts::VectorHelpers::phi;
0021
0022 namespace {
0023
0024 std::string effName(std::string_view label, std::string_view suffix) {
0025   return std::format("{}eff_{}", label, suffix);
0026 }
0027
0028 }  // namespace
0029
0030 namespace ActsExamples {
0031
0032 EffPlotTool::EffPlotTool(const EffPlotTool::Config& cfg,
0033                          Acts::Logging::Level lvl)
0034     : m_cfg(cfg), m_logger(Acts::getDefaultLogger("EffPlotTool", lvl)) {
0035   ACTS_DEBUG("Initialize the histograms for efficiency plots");
0036
0037   const std::string effTitlePtCut =
0038       std::format("{} efficiency with pT > {} GeV/c", m_cfg.label,
0039                   m_cfg.minTruthPt / Acts::UnitConstants::GeV);
0040
0041   // 1D efficiencies
0042   m_efficiencies1D.insert(
0043       {effName(m_cfg.label, "vs_eta"),
0044        Efficiency1(effName(m_cfg.label, "vs_eta"), effTitlePtCut,
0045                    std::array{m_cfg.varBinning.at("Eta")})});
0046   m_efficiencies1D.insert(
0047       {effName(m_cfg.label, "vs_phi"),
0048        Efficiency1(effName(m_cfg.label, "vs_phi"), effTitlePtCut,
0049                    std::array{m_cfg.varBinning.at("Phi")})});
0050   m_efficiencies1D.insert(
0051       {effName(m_cfg.label, "vs_pT"),
0052        Efficiency1(effName(m_cfg.label, "vs_pT"),
0053                    std::format("{} efficiency", m_cfg.label),
0054                    std::array{m_cfg.varBinning.at("Pt")})});
0055   m_efficiencies1D.insert(
0056       {effName(m_cfg.label, "vs_LogPt"),
0057        Efficiency1(effName(m_cfg.label, "vs_LogPt"),
0058                    std::format("{} efficiency", m_cfg.label),
0059                    std::array{m_cfg.varBinning.at("LogPt")})});
0060   m_efficiencies1D.insert(
0061       {effName(m_cfg.label, "vs_LowPt"),
0062        Efficiency1(effName(m_cfg.label, "vs_LowPt"),
0063                    std::format("{} efficiency", m_cfg.label),
0064                    std::array{m_cfg.varBinning.at("LowPt")})});
0065   m_efficiencies1D.insert(
0066       {effName(m_cfg.label, "vs_d0"),
0067        Efficiency1(effName(m_cfg.label, "vs_d0"), effTitlePtCut,
0068                    std::array{m_cfg.varBinning.at("D0")})});
0069   m_efficiencies1D.insert(
0070       {effName(m_cfg.label, "vs_z0"),
0071        Efficiency1(effName(m_cfg.label, "vs_z0"), effTitlePtCut,
0072                    std::array{m_cfg.varBinning.at("Z0")})});
0073   m_efficiencies1D.insert(
0074       {effName(m_cfg.label, "vs_DeltaR"),
0075        Efficiency1(effName(m_cfg.label, "vs_DeltaR"), effTitlePtCut,
0076                    std::array{m_cfg.varBinning.at("DeltaR")})});
0077   m_efficiencies1D.insert(
0078       {effName(m_cfg.label, "vs_prodR"),
0079        Efficiency1(effName(m_cfg.label, "vs_prodR"), effTitlePtCut,
0080                    std::array{m_cfg.varBinning.at("prodR")})});
0081
0082   // 2D efficiencies
0083   m_efficiencies2D.insert(
0084       {effName(m_cfg.label, "vs_eta_phi"),
0085        Efficiency2(effName(m_cfg.label, "vs_eta_phi"), effTitlePtCut,
0086                    std::array{m_cfg.varBinning.at("Eta"),
0087                               m_cfg.varBinning.at("Phi")})});
0088   m_efficiencies2D.insert(
0089       {effName(m_cfg.label, "vs_eta_pt"),
0090        Efficiency2(
0091            effName(m_cfg.label, "vs_eta_pt"),
0092            std::format("{} efficiency", m_cfg.label),
0093            std::array{m_cfg.varBinning.at("Eta"), m_cfg.varBinning.at("Pt")})});
0094
0095   const auto& etaAxis = m_cfg.varBinning.at("Eta");
0096   const auto& ptAxis = m_cfg.varBinning.at("Pt");
0097
0098   // efficiency vs eta in different pT ranges
0099   for (const auto& [i, ptRange] : Acts::enumerate(m_cfg.truthPtRangesForEta)) {
0100     const std::string name =
0101         std::format("{}eff_vs_eta_ptRange_{}", m_cfg.label, i);
0102     const std::string title =
0103         std::format("{} efficiency with pT in [{}, {}] GeV/c", m_cfg.label,
0104                     ptRange.first / Acts::UnitConstants::GeV,
0105                     ptRange.second / Acts::UnitConstants::GeV);
0106     m_trackEffVsEtaInPtRanges.emplace_back(name, title, std::array{etaAxis});
0107   }
0108
0109   // efficiency vs pT in different abs(eta) ranges
0110   for (const auto& [i, absEtaRange] :
0111        Acts::enumerate(m_cfg.truthAbsEtaRangesForPt)) {
0112     const std::string name =
0113         std::format("{}eff_vs_pT_absEtaRange_{}", m_cfg.label, i);
0114     const std::string title =
0115         std::format("{} efficiency with |#eta| in [{}, {}]", m_cfg.label,
0116                     absEtaRange.first, absEtaRange.second);
0117     m_trackEffVsPtInAbsEtaRanges.emplace_back(name, title, std::array{ptAxis});
0118   }
0119 }
0120
0121 void EffPlotTool::fill(const Acts::GeometryContext& gctx,
0122                        const SimParticleState& truthParticle,
0123                        const double deltaR, const bool status) {
0124   constexpr double nan = std::numeric_limits<double>::quiet_NaN();
0125
0126   const auto intersection =
0127       m_cfg.beamline
0128           ->intersect(gctx, truthParticle.position(), truthParticle.direction())
0129           .closest();
0130   Acts::Vector2 d0z0{nan, nan};
0131   if (intersection.isValid()) {
0132     auto localRes = m_cfg.beamline->globalToLocal(gctx, intersection.position(),
0133                                                   truthParticle.direction());
0134     if (localRes.ok()) {
0135       d0z0 = localRes.value();
0136     }
0137   }
0138
0139   const double t_phi = phi(truthParticle.direction());
0140   const double t_eta = eta(truthParticle.direction());
0141   const double t_absEta = std::abs(t_eta);
0142   const double t_pT = truthParticle.transverseMomentum();
0143   const double t_d0 = d0z0.x();
0144   const double t_z0 = d0z0.y();
0145   const double t_deltaR = deltaR;
0146   const double t_prodR = perp(truthParticle.position());
0147
0148   // cut on truth pT with the global range for the relevant plots
0149   if (t_pT >= m_cfg.minTruthPt) {
0150     m_efficiencies1D.at(effName(m_cfg.label, "vs_eta")).fill({t_eta}, status);
0151     m_efficiencies1D.at(effName(m_cfg.label, "vs_phi")).fill({t_phi}, status);
0152     m_efficiencies1D.at(effName(m_cfg.label, "vs_d0")).fill({t_d0}, status);
0153     m_efficiencies1D.at(effName(m_cfg.label, "vs_z0")).fill({t_z0}, status);
0154     m_efficiencies1D.at(effName(m_cfg.label, "vs_DeltaR"))
0155         .fill({t_deltaR}, status);
0156     m_efficiencies1D.at(effName(m_cfg.label, "vs_prodR"))
0157         .fill({t_prodR}, status);
0158
0159     m_efficiencies2D.at(effName(m_cfg.label, "vs_eta_phi"))
0160         .fill({t_eta, t_phi}, status);
0161   }
0162
0163   // do not cut on truth pT as it is a variable on the plot
0164   m_efficiencies1D.at(effName(m_cfg.label, "vs_pT")).fill({t_pT}, status);
0165   m_efficiencies1D.at(effName(m_cfg.label, "vs_LogPt")).fill({t_pT}, status);
0166   m_efficiencies1D.at(effName(m_cfg.label, "vs_LowPt")).fill({t_pT}, status);
0167   m_efficiencies2D.at(effName(m_cfg.label, "vs_eta_pt"))
0168       .fill({t_eta, t_pT}, status);
0169
0170   // fill the efficiency vs eta in different pT ranges
0171   for (auto&& [ptRange, eff] :
0172        Acts::zip(m_cfg.truthPtRangesForEta, m_trackEffVsEtaInPtRanges)) {
0173     if (t_pT >= ptRange.first && t_pT < ptRange.second) {
0174       eff.fill({t_eta}, status);
0175     }
0176   }
0177
0178   // fill the efficiency vs pT in different eta ranges
0179   for (auto&& [absEtaRange, eff] :
0180        Acts::zip(m_cfg.truthAbsEtaRangesForPt, m_trackEffVsPtInAbsEtaRanges)) {
0181     if (t_absEta >= absEtaRange.first && t_absEta < absEtaRange.second) {
0182       eff.fill({t_pT}, status);
0183     }
0184   }
0185 }
0186
0187 }  // namespace ActsExamples