EIC Software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
CKFPerformanceWriter.cpp
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file CKFPerformanceWriter.cpp
1 // This file is part of the Acts project.
2 //
3 // Copyright (C) 2020 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 
10 
15 
16 #include <numeric>
17 #include <stdexcept>
18 
19 #include <TFile.h>
20 #include <TTree.h>
21 
24  : WriterT(cfg.inputTrajectories, "CKFPerformanceWriter", lvl),
25  m_cfg(std::move(cfg)),
26  m_effPlotTool(m_cfg.effPlotToolConfig, lvl),
27  m_fakeRatePlotTool(m_cfg.fakeRatePlotToolConfig, lvl),
28  m_duplicationPlotTool(m_cfg.duplicationPlotToolConfig, lvl),
29  m_trackSummaryPlotTool(m_cfg.trackSummaryPlotToolConfig, lvl) {
30  // Input track and truth collection name
31  if (m_cfg.inputTrajectories.empty()) {
32  throw std::invalid_argument("Missing input trajectories collection");
33  }
34  if (m_cfg.inputParticles.empty()) {
35  throw std::invalid_argument("Missing input particles collection");
36  }
37  if (m_cfg.outputFilename.empty()) {
38  throw std::invalid_argument("Missing output filename");
39  }
40 
41  // the output file can not be given externally since TFile accesses to the
42  // same file from multiple threads are unsafe.
43  // must always be opened internally
45  m_outputFile = TFile::Open(path.c_str(), "RECREATE");
46  if (not m_outputFile) {
47  throw std::invalid_argument("Could not open '" + path + "'");
48  }
49 
50  // initialize the plot tools
55 }
56 
58  m_effPlotTool.clear(m_effPlotCache);
59  m_fakeRatePlotTool.clear(m_fakeRatePlotCache);
60  m_duplicationPlotTool.clear(m_duplicationPlotCache);
61  m_trackSummaryPlotTool.clear(m_trackSummaryPlotCache);
62  if (m_outputFile) {
63  m_outputFile->Close();
64  }
65 }
66 
68  if (m_outputFile) {
69  m_outputFile->cd();
70  m_effPlotTool.write(m_effPlotCache);
71  m_fakeRatePlotTool.write(m_fakeRatePlotCache);
72  m_duplicationPlotTool.write(m_duplicationPlotCache);
73  m_trackSummaryPlotTool.write(m_trackSummaryPlotCache);
74  ACTS_INFO("Wrote performance plots to '" << m_outputFile->GetPath() << "'");
75  }
76  return ProcessCode::SUCCESS;
77 }
78 
80  const AlgorithmContext& ctx, const TrajectoryContainer& trajectories) {
81  // The number of majority particle hits and fitted track parameters
82  using RecoTrackInfo = std::pair<size_t, Acts::BoundTrackParameters>;
83 
84  // Read truth particles from input collection
85  const auto& particles =
86  ctx.eventStore.get<SimParticleContainer>(m_cfg.inputParticles);
87 
88  // Exclusive access to the tree while writing
89  std::lock_guard<std::mutex> lock(m_writeMutex);
90 
91  // Counter of truth-matched reco tracks
92  std::map<ActsFatras::Barcode, std::vector<RecoTrackInfo>> matched;
93  // Counter of truth-unmatched reco tracks
94  std::map<ActsFatras::Barcode, size_t> unmatched;
95 
96  // Loop over all trajectories
97  for (const auto& traj : trajectories) {
98  // The trajectory entry indices and the multiTrajectory
99  const auto& [trackTips, mj] = traj.trajectory();
100  if (trackTips.empty()) {
101  ACTS_WARNING("Empty multiTrajectory.");
102  continue;
103  }
104 
105  // Loop over all trajectories in a multiTrajectory
106  for (const size_t& trackTip : trackTips) {
107  // Collect the trajectory summary info
108  auto trajState =
110  // Reco track selection
111  //@TODO: add interface for applying others cuts on reco tracks:
112  // -> pT, d0, z0, detector-specific hits/holes number cut
113  if (trajState.nMeasurements < m_cfg.nMeasurementsMin) {
114  continue;
115  }
116  // Check if the reco track has fitted track parameters
117  if (not traj.hasTrackParameters(trackTip)) {
118  ACTS_WARNING(
119  "No fitted track parameters for trajectory with entry index = "
120  << trackTip);
121  continue;
122  }
123  const auto& fittedParameters = traj.trackParameters(trackTip);
124  // Fill the trajectory summary info
125  m_trackSummaryPlotTool.fill(m_trackSummaryPlotCache, fittedParameters,
126  trajState.nStates, trajState.nMeasurements,
127  trajState.nOutliers, trajState.nHoles);
128 
129  // Get the majority truth particle to this track
130  std::vector<ParticleHitCount> particleHitCount =
131  traj.identifyMajorityParticle(trackTip);
132  if (particleHitCount.empty()) {
133  ACTS_WARNING(
134  "No truth particle associated with this trajectory with entry "
135  "index = "
136  << trackTip);
137  continue;
138  }
139  // Get the majority particleId and majority particle counts
140  // Note that the majority particle might be not in the truth seeds
141  // collection
142  ActsFatras::Barcode majorityParticleId =
143  particleHitCount.front().particleId;
144  size_t nMajorityHits = particleHitCount.front().hitCount;
145 
146  // Check if the trajectory is matched with truth.
147  // If not, it will be classified as 'fake'
148  bool isFake = false;
149  if (nMajorityHits * 1. / trajState.nMeasurements >=
150  m_cfg.truthMatchProbMin) {
151  matched[majorityParticleId].push_back(
152  {nMajorityHits, fittedParameters});
153  } else {
154  isFake = true;
155  unmatched[majorityParticleId]++;
156  }
157  // Fill fake rate plots
158  m_fakeRatePlotTool.fill(m_fakeRatePlotCache, fittedParameters, isFake);
159  } // end all trajectories in a multiTrajectory
160  } // end all multiTrajectories
161 
162  // Loop over all truth-matched reco tracks for duplication rate plots
163  for (auto& [particleId, matchedTracks] : matched) {
164  // Sort the reco tracks matched to this particle by the number of majority
165  // hits
166  std::sort(matchedTracks.begin(), matchedTracks.end(),
167  [](const RecoTrackInfo& lhs, const RecoTrackInfo& rhs) {
168  return lhs.first > rhs.first;
169  });
170  for (size_t itrack = 0; itrack < matchedTracks.size(); itrack++) {
171  const auto& [nMajorityHits, fittedParameters] = matchedTracks.at(itrack);
172  // The tracks with maximum number of majority hits is taken as the 'real'
173  // track; others are as 'duplicated'
174  bool isDuplicated = (itrack != 0);
175  // Fill the duplication rate
176  m_duplicationPlotTool.fill(m_duplicationPlotCache, fittedParameters,
177  isDuplicated);
178  }
179  }
180 
181  // Loop over all truth particle seeds for efficiency plots and reco details.
182  // These are filled w.r.t. truth particle seed info
183  for (const auto& particle : particles) {
184  auto particleId = particle.particleId();
185  // Investigate the truth-matched tracks
186  size_t nMatchedTracks = 0;
187  bool isReconstructed = false;
188  auto imatched = matched.find(particleId);
189  if (imatched != matched.end()) {
190  nMatchedTracks = imatched->second.size();
191  isReconstructed = true;
192  }
193  // Fill efficiency plots
194  m_effPlotTool.fill(m_effPlotCache, particle, isReconstructed);
195  // Fill number of duplicated tracks for this particle
196  m_duplicationPlotTool.fill(m_duplicationPlotCache, particle,
197  nMatchedTracks - 1);
198 
199  // Investigate the fake (i.e. truth-unmatched) tracks
200  size_t nFakeTracks = 0;
201  auto ifake = unmatched.find(particleId);
202  if (ifake != unmatched.end()) {
203  nFakeTracks = ifake->second;
204  }
205  // Fill number of reconstructed/truth-matched/fake tracks for this particle
206  m_fakeRatePlotTool.fill(m_fakeRatePlotCache, particle, nMatchedTracks,
207  nFakeTracks);
208  } // end all truth particles
209 
210  return ProcessCode::SUCCESS;
211 }