src/resolution/res_plot_tool.cpp

0001 /** TRACCC library, part of the ACTS project (R&D line)
0002  *
0003  * (c) 2022-2024 CERN for the benefit of the ACTS project
0004  *
0005  * Mozilla Public License Version 2.0
0006  */
0007
0008 // Library include(s).
0009 #include "res_plot_tool.hpp"
0010
0011 // Detray include(s).
0012 #include <detray/test/utils/statistics.hpp>
0013
0014 // ROOT include(s).
0015 #ifdef TRACCC_HAVE_ROOT
0016 #include <TF1.h>
0017 #include <TGraphErrors.h>
0018 #endif  // TRACCC_HAVE_ROOT
0019
0020 namespace traccc {
0021
0022 res_plot_tool::res_plot_tool(const res_plot_tool_config& cfg) : m_cfg(cfg) {}
0023
0024 void res_plot_tool::book(res_plot_cache& cache) const {
0025   plot_helpers::binning b_pull = m_cfg.var_binning.at("pull");
0026   plot_helpers::binning b_eta = m_cfg.var_binning.at("Eta");
0027   plot_helpers::binning b_pT = m_cfg.var_binning.at("Pt");
0028
0029   // Avoid unused variable warnings when building the code without ROOT.
0030   (void)cache;
0031   (void)b_pull;
0032   (void)b_eta;
0033   (void)b_pT;
0034
0035   for (std::size_t idx = 0; idx < m_cfg.param_names.size(); idx++) {
0036     std::string par_name = m_cfg.param_names.at(idx);
0037
0038     // Binning for residual is parameter dependent
0039     std::string par_residual = "residual_" + par_name;
0040     plot_helpers::binning b_residual = m_cfg.var_binning.at(par_residual);
0041
0042     // Avoid unused variable warnings when building the code without ROOT.
0043     (void)b_residual;
0044
0045 #ifdef TRACCC_HAVE_ROOT
0046     // residual distributions
0047     cache.residuals[par_name] = plot_helpers::book_histo(
0048         Form("res_%s", par_name.c_str()),
0049         Form("Residual of %s", par_name.c_str()), b_residual);
0050
0051     // pull distritutions
0052     cache.pulls[par_name] =
0053         plot_helpers::book_histo(Form("pull_%s", par_name.c_str()),
0054                                  Form("Pull of %s", par_name.c_str()), b_pull);
0055
0056     // resolution vs. eta plots
0057     cache.resolutions_eta[par_name] = plot_helpers::book_histo(
0058         Form("resolution_%s_vs_eta", par_name.c_str()),
0059         Form("Resolution of %s vs. eta", par_name.c_str()), b_eta);
0060
0061     // resolution vs. pT plots
0062     cache.resolutions_pT[par_name] = plot_helpers::book_histo(
0063         Form("resolution_%s_vs_pT", par_name.c_str()),
0064         Form("Resolution of %s vs. pT", par_name.c_str()), b_pT);
0065
0066     // residuals vs. eta plots
0067     cache.residuals_eta[par_name] = plot_helpers::book_histo(
0068         Form("residual_%s_vs_eta", par_name.c_str()),
0069         Form("Residual of %s vs. eta", par_name.c_str()), b_eta, b_residual);
0070
0071     // residuals vs. pT plots
0072     cache.residuals_pT[par_name] = plot_helpers::book_histo(
0073         Form("residual_%s_vs_pT", par_name.c_str()),
0074         Form("Residual of %s vs. pT", par_name.c_str()), b_pT, b_residual);
0075
0076     for (int i = 0; i < b_eta.n_bins; i++) {
0077       cache.residuals_per_eta[par_name][static_cast<unsigned long>(i)] =
0078           plot_helpers::book_histo(
0079               Form("residual_%s_at_%d_th_eta", par_name.c_str(), i),
0080               Form("Residual of %s at %d th eta", par_name.c_str(), i),
0081               b_residual);
0082     }
0083
0084     for (int i = 0; i < b_pT.n_bins; i++) {
0085       cache.residuals_per_pT[par_name][static_cast<unsigned long>(i)] =
0086           plot_helpers::book_histo(
0087               Form("residual_%s_at_%d_th_pT", par_name.c_str(), i),
0088               Form("Residual of %s at %d th pT", par_name.c_str(), i),
0089               b_residual);
0090     }
0091 #endif  // TRACCC_HAVE_ROOT
0092   }
0093 }
0094
0095 void res_plot_tool::fill(res_plot_cache& cache,
0096                          const bound_track_parameters<>& truth_param,
0097                          const bound_track_parameters<>& fit_param,
0098                          const particle& ptc) const {
0099   // Find index of eta and pT for resolution histogram
0100   const scalar eta = vector::eta(ptc.momentum);
0101   const scalar pT = std::hypot(ptc.momentum[0], ptc.momentum[1]);
0102
0103   // Avoid unused variable warnings when building the code without ROOT.
0104   (void)cache;
0105   (void)eta;
0106   (void)pT;
0107
0108   for (std::size_t idx = 0; idx < m_cfg.param_names.size(); idx++) {
0109     std::string par_name = m_cfg.param_names.at(idx);
0110
0111     scalar residual = 0.f;
0112     scalar pull = 0.f;
0113
0114     if (idx < e_bound_size) {
0115       residual = fit_param[idx] - truth_param[idx];
0116       pull = residual /
0117              std::sqrt(getter::element(fit_param.covariance(), idx, idx));
0118     } else if (par_name == "qopT") {
0119       residual = fit_param.qopT() - truth_param.qopT();
0120     } else if (par_name == "qopz") {
0121       residual = fit_param.qopz() - truth_param.qopz();
0122     }
0123
0124     // Avoid unused variable warnings when building the code without ROOT.
0125     (void)residual;
0126     (void)pull;
0127
0128 #ifdef TRACCC_HAVE_ROOT
0129     const auto eta_idx =
0130         std::clamp(cache.resolutions_eta[par_name]->FindBin(eta) - 1, 0,
0131                    cache.resolutions_eta[par_name]->GetNbinsX() - 1);
0132     const auto pT_idx =
0133         std::clamp(cache.resolutions_pT[par_name]->FindBin(pT) - 1, 0,
0134                    cache.resolutions_pT[par_name]->GetNbinsX() - 1);
0135
0136     cache.residuals.at(par_name)->Fill(residual);
0137     if (idx < e_bound_size) {
0138       cache.pulls.at(par_name)->Fill(pull);
0139     }
0140     cache.residuals_eta.at(par_name)->Fill(eta, residual);
0141     cache.residuals_pT.at(par_name)->Fill(pT, residual);
0142     cache.residuals_per_eta.at(par_name)
0143         .at(static_cast<std::size_t>(eta_idx))
0144         ->Fill(residual);
0145     cache.residuals_per_pT.at(par_name)
0146         .at(static_cast<std::size_t>(pT_idx))
0147         ->Fill(residual);
0148 #endif  // TRACCC_HAVE_ROOT
0149   }
0150 }
0151
0152 void res_plot_tool::write(res_plot_cache& cache) const {
0153   // Avoid unused variable warnings when building the code without ROOT.
0154   (void)cache;
0155
0156 #ifdef TRACCC_HAVE_ROOT
0157   for (const auto& residual : cache.residuals) {
0158     residual.second->Write();
0159   }
0160   for (const auto& pull : cache.pulls) {
0161     pull.second->Write();
0162   }
0163   for (const auto& residual : cache.residuals_eta) {
0164     residual.second->Write();
0165   }
0166   for (const auto& residual : cache.residuals_pT) {
0167     residual.second->Write();
0168   }
0169
0170   auto create_resolution_graph = [](std::shared_ptr<TH1> H, auto residuals,
0171                                     const char* x_axis_title,
0172                                     const char* y_axis_title) {
0173     const auto n_bins = H->GetNbinsX();
0174
0175     std::vector<float> sigmas;
0176
0177     for (int i = 0; i < n_bins; i++) {
0178       auto& data = residuals[static_cast<std::size_t>(i)];
0179
0180       // When there is no entry
0181       if (data->GetEntries() == 0) {
0182         H->SetBinContent(i, 0);
0183         sigmas.push_back(0);
0184         continue;
0185       }
0186
0187       // Function used for the fit.
0188       TF1 gaus{"gaus", "gaus", -1., 1.};
0189       double fit_par[3];
0190       auto res = data->Fit("gaus", "Q0S");
0191       gaus.GetParameters(&fit_par[0]);
0192       H->SetBinContent(i + 1, fit_par[2]);
0193       sigmas.push_back(static_cast<float>(gaus.GetParError(2)));
0194     }
0195
0196     std::unique_ptr<TGraphErrors> G = std::make_unique<TGraphErrors>(H.get());
0197
0198     for (int i = 0; i < n_bins; i++) {
0199       G->SetPointError(i, H->GetBinWidth(i) / 2.f,
0200                        sigmas[static_cast<std::size_t>(i)]);
0201     }
0202     G->SetName(H->GetName());
0203     G->SetMinimum(1e-5);
0204     G->SetMaximum(1e-1);
0205     G->GetHistogram()->GetYaxis()->SetMaxDigits(2);
0206     G->GetHistogram()->GetXaxis()->SetTitle(x_axis_title);
0207     G->GetHistogram()->GetYaxis()->SetTitle(y_axis_title);
0208
0209     return G;
0210   };
0211
0212   for (const auto& resolution : cache.resolutions_eta) {
0213     const auto& par_name = resolution.first;
0214     auto& H = resolution.second;
0215
0216     plot_helpers::binning binning =
0217         m_cfg.var_binning.at("resolution_" + par_name);
0218
0219     auto G = create_resolution_graph(H, cache.residuals_per_eta.at(par_name),
0220                                      "#eta", binning.title.c_str());
0221
0222     G->Write();
0223   }
0224
0225   for (const auto& resolution : cache.resolutions_pT) {
0226     const auto& par_name = resolution.first;
0227     auto& H = resolution.second;
0228
0229     plot_helpers::binning binning =
0230         m_cfg.var_binning.at("resolution_" + par_name);
0231
0232     auto G = create_resolution_graph(H, cache.residuals_per_pT.at(par_name),
0233                                      "p_{T} [GeV/c]", binning.title.c_str());
0234
0235     G->Write();
0236   }
0237
0238 #endif  // TRACCC_HAVE_ROOT
0239 }
0240
0241 }  // namespace traccc