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 "ActsExamples/Validation/ResPlotTool.hpp"
10 
11 #include "Acts/Surfaces/PerigeeSurface.hpp"
12 #include "Acts/Utilities/Helpers.hpp"
13 
14 ActsExamples::ResPlotTool::ResPlotTool(
15  const ActsExamples::ResPlotTool::Config& cfg, Acts::Logging::Level lvl)
16  : m_cfg(cfg), m_logger(Acts::getDefaultLogger("ResPlotTool", lvl)) {}
17 
18 void ActsExamples::ResPlotTool::book(
19  ResPlotTool::ResPlotCache& resPlotCache) const {
20  PlotHelpers::Binning bEta = m_cfg.varBinning.at("Eta");
21  PlotHelpers::Binning bPt = m_cfg.varBinning.at("Pt");
22  PlotHelpers::Binning bPull = m_cfg.varBinning.at("Pull");
23 
24  ACTS_DEBUG("Initialize the histograms for residual and pull plots");
25  for (unsigned int parID = 0; parID < Acts::eBoundSize; parID++) {
26  std::string parName = m_cfg.paramNames.at(parID);
27 
28  std::string parResidual = "Residual_" + parName;
29  // Binning for residual is parameter dependent
30  PlotHelpers::Binning bResidual = m_cfg.varBinning.at(parResidual);
31 
32  // residual distributions
33  resPlotCache.res[parName] = PlotHelpers::bookHisto(
34  Form("res_%s", parName.c_str()),
35  Form("Residual of %s", parName.c_str()), bResidual);
36  // residual vs eta scatter plots
37  resPlotCache.res_vs_eta[parName] = PlotHelpers::bookHisto(
38  Form("res_%s_vs_eta", parName.c_str()),
39  Form("Residual of %s vs eta", parName.c_str()), bEta, bResidual);
40  // residual mean in each eta bin
41  resPlotCache.resMean_vs_eta[parName] = PlotHelpers::bookHisto(
42  Form("resmean_%s_vs_eta", parName.c_str()),
43  Form("Residual mean of %s", parName.c_str()), bEta);
44  // residual width in each eta bin
45  resPlotCache.resWidth_vs_eta[parName] = PlotHelpers::bookHisto(
46  Form("reswidth_%s_vs_eta", parName.c_str()),
47  Form("Residual width of %s", parName.c_str()), bEta);
48  // residual vs pT scatter plots
49  resPlotCache.res_vs_pT[parName] = PlotHelpers::bookHisto(
50  Form("res_%s_vs_pT", parName.c_str()),
51  Form("Residual of %s vs pT", parName.c_str()), bPt, bResidual);
52  // residual mean in each pT bin
53  resPlotCache.resMean_vs_pT[parName] = PlotHelpers::bookHisto(
54  Form("resmean_%s_vs_pT", parName.c_str()),
55  Form("Residual mean of %s", parName.c_str()), bPt);
56  // residual width in each pT bin
57  resPlotCache.resWidth_vs_pT[parName] = PlotHelpers::bookHisto(
58  Form("reswidth_%s_vs_pT", parName.c_str()),
59  Form("Residual width of %s", parName.c_str()), bPt);
60 
61  // pull distritutions
62  resPlotCache.pull[parName] =
63  PlotHelpers::bookHisto(Form("pull_%s", parName.c_str()),
64  Form("Pull of %s", parName.c_str()), bPull);
65  // pull vs eta scatter plots
66  resPlotCache.pull_vs_eta[parName] = PlotHelpers::bookHisto(
67  Form("pull_%s_vs_eta", parName.c_str()),
68  Form("Pull of %s vs eta", parName.c_str()), bEta, bPull);
69  // pull mean in each eta bin
70  resPlotCache.pullMean_vs_eta[parName] =
71  PlotHelpers::bookHisto(Form("pullmean_%s_vs_eta", parName.c_str()),
72  Form("Pull mean of %s", parName.c_str()), bEta);
73  // pull width in each eta bin
74  resPlotCache.pullWidth_vs_eta[parName] =
75  PlotHelpers::bookHisto(Form("pullwidth_%s_vs_eta", parName.c_str()),
76  Form("Pull width of %s", parName.c_str()), bEta);
77  // pull vs pT scatter plots
78  resPlotCache.pull_vs_pT[parName] = PlotHelpers::bookHisto(
79  Form("pull_%s_vs_pT", parName.c_str()),
80  Form("Pull of %s vs pT", parName.c_str()), bPt, bPull);
81  // pull mean in each pT bin
82  resPlotCache.pullMean_vs_pT[parName] =
83  PlotHelpers::bookHisto(Form("pullmean_%s_vs_pT", parName.c_str()),
84  Form("Pull mean of %s", parName.c_str()), bPt);
85  // pull width in each pT bin
86  resPlotCache.pullWidth_vs_pT[parName] =
87  PlotHelpers::bookHisto(Form("pullwidth_%s_vs_pT", parName.c_str()),
88  Form("Pull width of %s", parName.c_str()), bPt);
89  }
90 }
91 
92 void ActsExamples::ResPlotTool::clear(ResPlotCache& resPlotCache) const {
93  ACTS_DEBUG("Delete the hists.");
94  for (unsigned int parID = 0; parID < Acts::eBoundSize; parID++) {
95  std::string parName = m_cfg.paramNames.at(parID);
96  delete resPlotCache.res.at(parName);
97  delete resPlotCache.res_vs_eta.at(parName);
98  delete resPlotCache.resMean_vs_eta.at(parName);
99  delete resPlotCache.resWidth_vs_eta.at(parName);
100  delete resPlotCache.res_vs_pT.at(parName);
101  delete resPlotCache.resMean_vs_pT.at(parName);
102  delete resPlotCache.resWidth_vs_pT.at(parName);
103  delete resPlotCache.pull.at(parName);
104  delete resPlotCache.pull_vs_eta.at(parName);
105  delete resPlotCache.pullMean_vs_eta.at(parName);
106  delete resPlotCache.pullWidth_vs_eta.at(parName);
107  delete resPlotCache.pull_vs_pT.at(parName);
108  delete resPlotCache.pullMean_vs_pT.at(parName);
109  delete resPlotCache.pullWidth_vs_pT.at(parName);
110  }
111 }
112 
113 void ActsExamples::ResPlotTool::write(
114  const ResPlotTool::ResPlotCache& resPlotCache) const {
115  ACTS_DEBUG("Write the hists to output file.");
116  for (unsigned int parID = 0; parID < Acts::eBoundSize; parID++) {
117  std::string parName = m_cfg.paramNames.at(parID);
118  resPlotCache.res.at(parName)->Write();
119  resPlotCache.res_vs_eta.at(parName)->Write();
120  resPlotCache.resMean_vs_eta.at(parName)->Write();
121  resPlotCache.resWidth_vs_eta.at(parName)->Write();
122  resPlotCache.res_vs_pT.at(parName)->Write();
123  resPlotCache.resMean_vs_pT.at(parName)->Write();
124  resPlotCache.resWidth_vs_pT.at(parName)->Write();
125  resPlotCache.pull.at(parName)->Write();
126  resPlotCache.pull_vs_eta.at(parName)->Write();
127  resPlotCache.pullMean_vs_eta.at(parName)->Write();
128  resPlotCache.pullWidth_vs_eta.at(parName)->Write();
129  resPlotCache.pull_vs_pT.at(parName)->Write();
130  resPlotCache.pullMean_vs_pT.at(parName)->Write();
131  resPlotCache.pullWidth_vs_pT.at(parName)->Write();
132  }
133 }
134 
135 void ActsExamples::ResPlotTool::fill(
136  ResPlotTool::ResPlotCache& resPlotCache, const Acts::GeometryContext& gctx,
137  const ActsFatras::Particle& truthParticle,
138  const Acts::BoundTrackParameters& fittedParamters) const {
139  using ParametersVector = Acts::BoundTrackParameters::ParametersVector;
140  using Acts::VectorHelpers::eta;
141  using Acts::VectorHelpers::perp;
142  using Acts::VectorHelpers::phi;
143  using Acts::VectorHelpers::theta;
144 
145  // get the fitted parameter (at perigee surface) and its error
146  auto trackParameter = fittedParamters.parameters();
147  auto covariance = *fittedParamters.covariance();
148 
149  // get the perigee surface
150  auto pSurface = &fittedParamters.referenceSurface();
151 
152  // get the truth position and momentum
153  ParametersVector truthParameter = ParametersVector::Zero();
154 
155  // get the truth perigee parameter
156  auto lpResult = pSurface->globalToLocal(gctx, truthParticle.position(),
157  truthParticle.unitDirection());
158  if (lpResult.ok()) {
159  truthParameter[Acts::BoundIndices::eBoundLoc0] =
160  lpResult.value()[Acts::BoundIndices::eBoundLoc0];
161  truthParameter[Acts::BoundIndices::eBoundLoc1] =
162  lpResult.value()[Acts::BoundIndices::eBoundLoc1];
163  } else {
164  ACTS_ERROR("Global to local transformation did not succeed.");
165  }
166  truthParameter[Acts::BoundIndices::eBoundPhi] =
167  phi(truthParticle.unitDirection());
168  truthParameter[Acts::BoundIndices::eBoundTheta] =
169  theta(truthParticle.unitDirection());
170  truthParameter[Acts::BoundIndices::eBoundQOverP] =
171  truthParticle.charge() / truthParticle.absMomentum();
172  truthParameter[Acts::BoundIndices::eBoundTime] = truthParticle.time();
173 
174  // get the truth eta and pT
175  const auto truthEta = eta(truthParticle.unitDirection());
176  const auto truthPt = truthParticle.transverseMomentum();
177 
178  // fill the histograms for residual and pull
179  for (unsigned int parID = 0; parID < Acts::eBoundSize; parID++) {
180  std::string parName = m_cfg.paramNames.at(parID);
181  float residual = trackParameter[parID] - truthParameter[parID];
182  PlotHelpers::fillHisto(resPlotCache.res.at(parName), residual);
183  PlotHelpers::fillHisto(resPlotCache.res_vs_eta.at(parName), truthEta,
184  residual);
185  PlotHelpers::fillHisto(resPlotCache.res_vs_pT.at(parName), truthPt,
186  residual);
187  if (covariance(parID, parID) > 0) {
188  float pull = residual / sqrt(covariance(parID, parID));
189  PlotHelpers::fillHisto(resPlotCache.pull[parName], pull);
190  PlotHelpers::fillHisto(resPlotCache.pull_vs_eta.at(parName), truthEta,
191  pull);
192  PlotHelpers::fillHisto(resPlotCache.pull_vs_pT.at(parName), truthPt,
193  pull);
194  } else {
195  ACTS_WARNING("Fitted track parameter :" << parName
196  << " has negative covariance = "
197  << covariance(parID, parID));
198  }
199  }
200 }
201 
202 // get the mean and width of residual/pull in each eta/pT bin and fill them into
203 // histograms
204 void ActsExamples::ResPlotTool::refinement(
205  ResPlotTool::ResPlotCache& resPlotCache) const {
206  PlotHelpers::Binning bEta = m_cfg.varBinning.at("Eta");
207  PlotHelpers::Binning bPt = m_cfg.varBinning.at("Pt");
208  for (unsigned int parID = 0; parID < Acts::eBoundSize; parID++) {
209  std::string parName = m_cfg.paramNames.at(parID);
210  // refine the plots vs eta
211  for (int j = 1; j <= bEta.nBins; j++) {
212  TH1D* temp_res = resPlotCache.res_vs_eta.at(parName)->ProjectionY(
213  Form("%s_projy_bin%d", "Residual_vs_eta_Histo", j), j, j);
214  PlotHelpers::anaHisto(temp_res, j,
215  resPlotCache.resMean_vs_eta.at(parName),
216  resPlotCache.resWidth_vs_eta.at(parName));
217 
218  TH1D* temp_pull = resPlotCache.pull_vs_eta.at(parName)->ProjectionY(
219  Form("%s_projy_bin%d", "Pull_vs_eta_Histo", j), j, j);
220  PlotHelpers::anaHisto(temp_pull, j,
221  resPlotCache.pullMean_vs_eta.at(parName),
222  resPlotCache.pullWidth_vs_eta.at(parName));
223  }
224 
225  // refine the plots vs pT
226  for (int j = 1; j <= bPt.nBins; j++) {
227  TH1D* temp_res = resPlotCache.res_vs_pT.at(parName)->ProjectionY(
228  Form("%s_projy_bin%d", "Residual_vs_pT_Histo", j), j, j);
229  PlotHelpers::anaHisto(temp_res, j, resPlotCache.resMean_vs_pT.at(parName),
230  resPlotCache.resWidth_vs_pT.at(parName));
231 
232  TH1D* temp_pull = resPlotCache.pull_vs_pT.at(parName)->ProjectionY(
233  Form("%s_projy_bin%d", "Pull_vs_pT_Histo", j), j, j);
234  PlotHelpers::anaHisto(temp_pull, j,
235  resPlotCache.pullMean_vs_pT.at(parName),
236  resPlotCache.pullWidth_vs_pT.at(parName));
237  }
238  }
239 }