ECCE @ EIC Software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
ResPlotTool.cpp
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file ResPlotTool.cpp
1 // This file is part of the Acts project.
2 //
3 // Copyright (C) 2019 CERN for the benefit of the Acts project
4 //
5 // This Source Code Form is subject to the terms of the Mozilla Public
6 // License, v. 2.0. If a copy of the MPL was not distributed with this
7 // file, You can obtain one at http://mozilla.org/MPL/2.0/.
8 
9 #include "ACTFW/Validation/ResPlotTool.hpp"
10 
11 #include "Acts/Surfaces/PerigeeSurface.hpp"
12 #include "Acts/Utilities/Helpers.hpp"
13 
14 FW::ResPlotTool::ResPlotTool(const FW::ResPlotTool::Config& cfg,
15  Acts::Logging::Level lvl)
16  : m_cfg(cfg), m_logger(Acts::getDefaultLogger("ResPlotTool", lvl)) {}
17 
18 void FW::ResPlotTool::book(ResPlotTool::ResPlotCache& resPlotCache) const {
19  PlotHelpers::Binning bEta = m_cfg.varBinning.at("Eta");
20  PlotHelpers::Binning bPt = m_cfg.varBinning.at("Pt");
21  PlotHelpers::Binning bPull = m_cfg.varBinning.at("Pull");
22 
23  ACTS_DEBUG("Initialize the histograms for residual and pull plots");
24  for (unsigned int parID = 0; parID < Acts::eBoundParametersSize; parID++) {
25  std::string parName = m_cfg.paramNames.at(parID);
26 
27  std::string parResidual = "Residual_" + parName;
28  // Binning for residual is parameter dependent
29  PlotHelpers::Binning bResidual = m_cfg.varBinning.at(parResidual);
30 
31  // residual distributions
32  resPlotCache.res[parName] = PlotHelpers::bookHisto(
33  Form("res_%s", parName.c_str()),
34  Form("Residual of %s", parName.c_str()), bResidual);
35  // residual vs eta scatter plots
36  resPlotCache.res_vs_eta[parName] = PlotHelpers::bookHisto(
37  Form("res_%s_vs_eta", parName.c_str()),
38  Form("Residual of %s vs eta", parName.c_str()), bEta, bResidual);
39  // residual mean in each eta bin
40  resPlotCache.resmean_vs_eta[parName] = PlotHelpers::bookHisto(
41  Form("resmean_%s_vs_eta", parName.c_str()),
42  Form("Residual mean of %s", parName.c_str()), bEta);
43  // residual width in each eta bin
44  resPlotCache.reswidth_vs_eta[parName] = PlotHelpers::bookHisto(
45  Form("reswidth_%s_vs_eta", parName.c_str()),
46  Form("Residual width of %s", parName.c_str()), bEta);
47  // residual vs pT scatter plots
48  resPlotCache.res_vs_pT[parName] = PlotHelpers::bookHisto(
49  Form("res_%s_vs_pT", parName.c_str()),
50  Form("Residual of %s vs pT", parName.c_str()), bPt, bResidual);
51  // residual mean in each pT bin
52  resPlotCache.resmean_vs_pT[parName] = PlotHelpers::bookHisto(
53  Form("resmean_%s_vs_pT", parName.c_str()),
54  Form("Residual mean of %s", parName.c_str()), bPt);
55  // residual width in each pT bin
56  resPlotCache.reswidth_vs_pT[parName] = PlotHelpers::bookHisto(
57  Form("reswidth_%s_vs_pT", parName.c_str()),
58  Form("Residual width of %s", parName.c_str()), bPt);
59 
60  // pull distritutions
61  resPlotCache.pull[parName] =
62  PlotHelpers::bookHisto(Form("pull_%s", parName.c_str()),
63  Form("Pull of %s", parName.c_str()), bPull);
64  // pull vs eta scatter plots
65  resPlotCache.pull_vs_eta[parName] = PlotHelpers::bookHisto(
66  Form("pull_%s_vs_eta", parName.c_str()),
67  Form("Pull of %s vs eta", parName.c_str()), bEta, bPull);
68  // pull mean in each eta bin
69  resPlotCache.pullmean_vs_eta[parName] =
70  PlotHelpers::bookHisto(Form("pullmean_%s_vs_eta", parName.c_str()),
71  Form("Pull mean of %s", parName.c_str()), bEta);
72  // pull width in each eta bin
73  resPlotCache.pullwidth_vs_eta[parName] =
74  PlotHelpers::bookHisto(Form("pullwidth_%s_vs_eta", parName.c_str()),
75  Form("Pull width of %s", parName.c_str()), bEta);
76  // pull vs pT scatter plots
77  resPlotCache.pull_vs_pT[parName] = PlotHelpers::bookHisto(
78  Form("pull_%s_vs_pT", parName.c_str()),
79  Form("Pull of %s vs pT", parName.c_str()), bPt, bPull);
80  // pull mean in each pT bin
81  resPlotCache.pullmean_vs_pT[parName] =
82  PlotHelpers::bookHisto(Form("pullmean_%s_vs_pT", parName.c_str()),
83  Form("Pull mean of %s", parName.c_str()), bPt);
84  // pull width in each pT bin
85  resPlotCache.pullwidth_vs_pT[parName] =
86  PlotHelpers::bookHisto(Form("pullwidth_%s_vs_pT", parName.c_str()),
87  Form("Pull width of %s", parName.c_str()), bPt);
88  }
89 }
90 
91 void FW::ResPlotTool::clear(ResPlotCache& resPlotCache) const {
92  ACTS_DEBUG("Delete the hists.");
93  for (unsigned int parID = 0; parID < Acts::eBoundParametersSize; parID++) {
94  std::string parName = m_cfg.paramNames.at(parID);
95  delete resPlotCache.res.at(parName);
96  delete resPlotCache.res_vs_eta.at(parName);
97  delete resPlotCache.resmean_vs_eta.at(parName);
98  delete resPlotCache.reswidth_vs_eta.at(parName);
99  delete resPlotCache.res_vs_pT.at(parName);
100  delete resPlotCache.resmean_vs_pT.at(parName);
101  delete resPlotCache.reswidth_vs_pT.at(parName);
102  delete resPlotCache.pull.at(parName);
103  delete resPlotCache.pull_vs_eta.at(parName);
104  delete resPlotCache.pullmean_vs_eta.at(parName);
105  delete resPlotCache.pullwidth_vs_eta.at(parName);
106  delete resPlotCache.pull_vs_pT.at(parName);
107  delete resPlotCache.pullmean_vs_pT.at(parName);
108  delete resPlotCache.pullwidth_vs_pT.at(parName);
109  }
110 }
111 
112 void FW::ResPlotTool::write(
113  const ResPlotTool::ResPlotCache& resPlotCache) const {
114  ACTS_DEBUG("Write the hists to output file.");
115  for (unsigned int parID = 0; parID < Acts::eBoundParametersSize; parID++) {
116  std::string parName = m_cfg.paramNames.at(parID);
117  resPlotCache.res.at(parName)->Write();
118  resPlotCache.res_vs_eta.at(parName)->Write();
119  resPlotCache.resmean_vs_eta.at(parName)->Write();
120  resPlotCache.reswidth_vs_eta.at(parName)->Write();
121  resPlotCache.res_vs_pT.at(parName)->Write();
122  resPlotCache.resmean_vs_pT.at(parName)->Write();
123  resPlotCache.reswidth_vs_pT.at(parName)->Write();
124  resPlotCache.pull.at(parName)->Write();
125  resPlotCache.pull_vs_eta.at(parName)->Write();
126  resPlotCache.pullmean_vs_eta.at(parName)->Write();
127  resPlotCache.pullwidth_vs_eta.at(parName)->Write();
128  resPlotCache.pull_vs_pT.at(parName)->Write();
129  resPlotCache.pullmean_vs_pT.at(parName)->Write();
130  resPlotCache.pullwidth_vs_pT.at(parName)->Write();
131  }
132 }
133 
134 void FW::ResPlotTool::fill(ResPlotTool::ResPlotCache& resPlotCache,
135  const Acts::GeometryContext& gctx,
136  const ActsFatras::Particle& truthParticle,
137  const Acts::BoundParameters& fittedParamters) const {
138  using Acts::VectorHelpers::eta;
139  using Acts::VectorHelpers::perp;
140  using Acts::VectorHelpers::phi;
141  using Acts::VectorHelpers::theta;
142 
143  // get the fitted parameter (at perigee surface) and its error
144  auto trackParameter = fittedParamters.parameters();
145  auto covariance = *fittedParamters.covariance();
146 
147  // get the perigee surface
148  auto pSurface = &fittedParamters.referenceSurface();
149 
150  // get the truth position and momentum
151  ParVector_t truthParameter;
152 
153  // get the truth perigee parameter
154  Acts::Vector2D local(0., 0.);
155  pSurface->globalToLocal(gctx, truthParticle.position(),
156  truthParticle.unitDirection(), local);
157  truthParameter[Acts::ParDef::eLOC_D0] = local[Acts::ParDef::eLOC_D0];
158  truthParameter[Acts::ParDef::eLOC_Z0] = local[Acts::ParDef::eLOC_Z0];
159  truthParameter[Acts::ParDef::ePHI] = phi(truthParticle.unitDirection());
160  truthParameter[Acts::ParDef::eTHETA] = theta(truthParticle.unitDirection());
161  truthParameter[Acts::ParDef::eQOP] =
162  truthParticle.charge() / truthParticle.absMomentum();
163  truthParameter[Acts::ParDef::eT] = truthParticle.time();
164 
165  // get the truth eta and pT
166  const auto truthEta = eta(truthParticle.unitDirection());
167  const auto truthPt = truthParticle.transverseMomentum();
168 
169  // fill the histograms for residual and pull
170  for (unsigned int parID = 0; parID < Acts::eBoundParametersSize; parID++) {
171  std::string parName = m_cfg.paramNames.at(parID);
172  float residual = trackParameter[parID] - truthParameter[parID];
173  PlotHelpers::fillHisto(resPlotCache.res.at(parName), residual);
174  PlotHelpers::fillHisto(resPlotCache.res_vs_eta.at(parName), truthEta,
175  residual);
176  PlotHelpers::fillHisto(resPlotCache.res_vs_pT.at(parName), truthPt,
177  residual);
178  if (covariance(parID, parID) > 0) {
179  float pull = residual / sqrt(covariance(parID, parID));
180  PlotHelpers::fillHisto(resPlotCache.pull[parName], pull);
181  PlotHelpers::fillHisto(resPlotCache.pull_vs_eta.at(parName), truthEta,
182  pull);
183  PlotHelpers::fillHisto(resPlotCache.pull_vs_pT.at(parName), truthPt,
184  pull);
185  } else {
186  ACTS_WARNING("Fitted track parameter :" << parName
187  << " has negative covariance = "
188  << covariance(parID, parID));
189  }
190  }
191 }
192 
193 // get the mean and width of residual/pull in each eta/pT bin and fill them into
194 // histograms
195 void FW::ResPlotTool::refinement(
196  ResPlotTool::ResPlotCache& resPlotCache) const {
197  PlotHelpers::Binning bEta = m_cfg.varBinning.at("Eta");
198  PlotHelpers::Binning bPt = m_cfg.varBinning.at("Pt");
199  for (unsigned int parID = 0; parID < Acts::eBoundParametersSize; parID++) {
200  std::string parName = m_cfg.paramNames.at(parID);
201  // refine the plots vs eta
202  for (int j = 1; j <= bEta.nBins; j++) {
203  TH1D* temp_res = resPlotCache.res_vs_eta.at(parName)->ProjectionY(
204  Form("%s_projy_bin%d", "Residual_vs_eta_Histo", j), j, j);
205  PlotHelpers::anaHisto(temp_res, j,
206  resPlotCache.resmean_vs_eta.at(parName),
207  resPlotCache.reswidth_vs_eta.at(parName));
208 
209  TH1D* temp_pull = resPlotCache.pull_vs_eta.at(parName)->ProjectionY(
210  Form("%s_projy_bin%d", "Pull_vs_eta_Histo", j), j, j);
211  PlotHelpers::anaHisto(temp_pull, j,
212  resPlotCache.pullmean_vs_eta.at(parName),
213  resPlotCache.pullwidth_vs_eta.at(parName));
214  }
215 
216  // refine the plots vs pT
217  for (int j = 1; j <= bPt.nBins; j++) {
218  TH1D* temp_res = resPlotCache.res_vs_pT.at(parName)->ProjectionY(
219  Form("%s_projy_bin%d", "Residual_vs_pT_Histo", j), j, j);
220  PlotHelpers::anaHisto(temp_res, j, resPlotCache.resmean_vs_pT.at(parName),
221  resPlotCache.reswidth_vs_pT.at(parName));
222 
223  TH1D* temp_pull = resPlotCache.pull_vs_pT.at(parName)->ProjectionY(
224  Form("%s_projy_bin%d", "Pull_vs_pT_Histo", j), j, j);
225  PlotHelpers::anaHisto(temp_pull, j,
226  resPlotCache.pullmean_vs_pT.at(parName),
227  resPlotCache.pullwidth_vs_pT.at(parName));
228  }
229  }
230 }