Analysis Software
Documentation for sPHENIX simulation software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
RootTrajectorySummaryReader.cpp
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file RootTrajectorySummaryReader.cpp
1 // This file is part of the Acts project.
2 //
3 // Copyright (C) 2021 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/Io/Root/RootTrajectorySummaryReader.hpp"
10 
11 #include "Acts/Definitions/Algebra.hpp"
12 #include "Acts/Definitions/TrackParametrization.hpp"
13 #include "Acts/EventData/TrackParameters.hpp"
14 #include "Acts/Surfaces/PerigeeSurface.hpp"
15 #include "Acts/Surfaces/Surface.hpp"
16 #include "Acts/Utilities/Logger.hpp"
17 #include "ActsExamples/EventData/SimParticle.hpp"
18 #include "ActsExamples/Framework/AlgorithmContext.hpp"
19 #include "ActsFatras/EventData/Particle.hpp"
20 
21 #include <iostream>
22 #include <stdexcept>
23 
24 #include <TChain.h>
25 #include <TMathBase.h>
26 
27 ActsExamples::RootTrajectorySummaryReader::RootTrajectorySummaryReader(
28  const ActsExamples::RootTrajectorySummaryReader::Config& config,
29  Acts::Logging::Level level)
30  : ActsExamples::IReader(),
31  m_logger{Acts::getDefaultLogger(name(), level)},
32  m_cfg(config) {
33  m_inputChain = new TChain(m_cfg.treeName.c_str());
34 
35  if (m_cfg.filePath.empty()) {
36  throw std::invalid_argument("Missing input filename");
37  }
38 
39  m_outputTrackParameters.initialize(m_cfg.outputTracks);
40  m_outputParticles.initialize(m_cfg.outputParticles);
41 
42  // Set the branches
43  m_inputChain->SetBranchAddress("event_nr", &m_eventNr);
44  m_inputChain->SetBranchAddress("multiTraj_nr", &m_multiTrajNr);
45  m_inputChain->SetBranchAddress("subTraj_nr", &m_subTrajNr);
46 
47  // These info is not really stored in the event store, but still read in
48  m_inputChain->SetBranchAddress("nStates", &m_nStates);
49  m_inputChain->SetBranchAddress("nMeasurements", &m_nMeasurements);
50  m_inputChain->SetBranchAddress("nOutliers", &m_nOutliers);
51  m_inputChain->SetBranchAddress("nHoles", &m_nHoles);
52  m_inputChain->SetBranchAddress("chi2Sum", &m_chi2Sum);
53  m_inputChain->SetBranchAddress("NDF", &m_NDF);
54  m_inputChain->SetBranchAddress("measurementChi2", &m_measurementChi2);
55  m_inputChain->SetBranchAddress("outlierChi2", &m_outlierChi2);
56  m_inputChain->SetBranchAddress("measurementVolume", &m_measurementVolume);
57  m_inputChain->SetBranchAddress("measurementLayer", &m_measurementLayer);
58  m_inputChain->SetBranchAddress("outlierVolume", &m_outlierVolume);
59  m_inputChain->SetBranchAddress("outlierLayer", &m_outlierLayer);
60 
61  m_inputChain->SetBranchAddress("majorityParticleId", &m_majorityParticleId);
62  m_inputChain->SetBranchAddress("nMajorityHits", &m_nMajorityHits);
63  m_inputChain->SetBranchAddress("t_charge", &m_t_charge);
64  m_inputChain->SetBranchAddress("t_time", &m_t_time);
65  m_inputChain->SetBranchAddress("t_vx", &m_t_vx);
66  m_inputChain->SetBranchAddress("t_vy", &m_t_vy);
67  m_inputChain->SetBranchAddress("t_vz", &m_t_vz);
68  m_inputChain->SetBranchAddress("t_px", &m_t_px);
69  m_inputChain->SetBranchAddress("t_py", &m_t_py);
70  m_inputChain->SetBranchAddress("t_pz", &m_t_pz);
71  m_inputChain->SetBranchAddress("t_theta", &m_t_theta);
72  m_inputChain->SetBranchAddress("t_phi", &m_t_phi);
73  m_inputChain->SetBranchAddress("t_eta", &m_t_eta);
74  m_inputChain->SetBranchAddress("t_pT", &m_t_pT);
75 
76  m_inputChain->SetBranchAddress("hasFittedParams", &m_hasFittedParams);
77  m_inputChain->SetBranchAddress("eLOC0_fit", &m_eLOC0_fit);
78  m_inputChain->SetBranchAddress("eLOC1_fit", &m_eLOC1_fit);
79  m_inputChain->SetBranchAddress("ePHI_fit", &m_ePHI_fit);
80  m_inputChain->SetBranchAddress("eTHETA_fit", &m_eTHETA_fit);
81  m_inputChain->SetBranchAddress("eQOP_fit", &m_eQOP_fit);
82  m_inputChain->SetBranchAddress("eT_fit", &m_eT_fit);
83  m_inputChain->SetBranchAddress("err_eLOC0_fit", &m_err_eLOC0_fit);
84  m_inputChain->SetBranchAddress("err_eLOC1_fit", &m_err_eLOC1_fit);
85  m_inputChain->SetBranchAddress("err_ePHI_fit", &m_err_ePHI_fit);
86  m_inputChain->SetBranchAddress("err_eTHETA_fit", &m_err_eTHETA_fit);
87  m_inputChain->SetBranchAddress("err_eQOP_fit", &m_err_eQOP_fit);
88  m_inputChain->SetBranchAddress("err_eT_fit", &m_err_eT_fit);
89 
90  auto path = m_cfg.filePath;
91 
92  // add file to the input chain
93  m_inputChain->Add(path.c_str());
94  ACTS_DEBUG("Adding File " << path << " to tree '" << m_cfg.treeName << "'.");
95 
96  m_events = m_inputChain->GetEntries();
97  ACTS_DEBUG("The full chain has " << m_events << " entries.");
98 
99  // If the events are not in order, get the entry numbers for ordered events
100  if (not m_cfg.orderedEvents) {
101  m_entryNumbers.resize(m_events);
102  m_inputChain->Draw("event_nr", "", "goff");
103  // Sort to get the entry numbers of the ordered events
104  TMath::Sort(m_inputChain->GetEntries(), m_inputChain->GetV1(),
105  m_entryNumbers.data(), false);
106  }
107 }
108 
109 std::pair<size_t, size_t>
110 ActsExamples::RootTrajectorySummaryReader::availableEvents() const {
111  return {0u, m_events};
112 }
113 
114 ActsExamples::RootTrajectorySummaryReader::~RootTrajectorySummaryReader() {
115  delete m_multiTrajNr;
116  delete m_subTrajNr;
117  delete m_nStates;
118  delete m_nMeasurements;
119  delete m_nOutliers;
120  delete m_nHoles;
121  delete m_chi2Sum;
122  delete m_NDF;
123  delete m_measurementChi2;
124  delete m_outlierChi2;
125  delete m_measurementVolume;
126  delete m_measurementLayer;
127  delete m_outlierVolume;
128  delete m_outlierLayer;
129  delete m_majorityParticleId;
130  delete m_nMajorityHits;
131  delete m_t_charge;
132  delete m_t_time;
133  delete m_t_vx;
134  delete m_t_vy;
135  delete m_t_vz;
136  delete m_t_px;
137  delete m_t_py;
138  delete m_t_pz;
139  delete m_t_theta;
140  delete m_t_phi;
141  delete m_t_pT;
142  delete m_t_eta;
143  delete m_hasFittedParams;
144  delete m_eLOC0_fit;
145  delete m_eLOC1_fit;
146  delete m_ePHI_fit;
147  delete m_eTHETA_fit;
148  delete m_eQOP_fit;
149  delete m_eT_fit;
150  delete m_err_eLOC0_fit;
151  delete m_err_eLOC1_fit;
152  delete m_err_ePHI_fit;
153  delete m_err_eTHETA_fit;
154  delete m_err_eQOP_fit;
155  delete m_err_eT_fit;
156 }
157 
158 ActsExamples::ProcessCode ActsExamples::RootTrajectorySummaryReader::read(
159  const ActsExamples::AlgorithmContext& context) {
160  ACTS_DEBUG("Trying to read recorded tracks.");
161 
162  // read in the fitted track parameters and particles
163  if (m_inputChain != nullptr && context.eventNumber < m_events) {
164  // lock the mutex
165  std::lock_guard<std::mutex> lock(m_read_mutex);
166  // now read
167 
168  std::shared_ptr<Acts::PerigeeSurface> perigeeSurface =
169  Acts::Surface::makeShared<Acts::PerigeeSurface>(
170  Acts::Vector3(0., 0., 0.));
171 
172  // The collection to be written
173  std::vector<Acts::BoundTrackParameters> trackParameterCollection;
174  SimParticleContainer truthParticleCollection;
175 
176  // Read the correct entry
177  auto entry = context.eventNumber;
178  if (not m_cfg.orderedEvents and entry < m_entryNumbers.size()) {
179  entry = m_entryNumbers[entry];
180  }
181  m_inputChain->GetEntry(entry);
182  ACTS_INFO("Reading event: " << context.eventNumber
183  << " stored as entry: " << entry);
184 
185  unsigned int nTracks = m_eLOC0_fit->size();
186  for (unsigned int i = 0; i < nTracks; i++) {
187  Acts::BoundVector paramVec;
188  paramVec << (*m_eLOC0_fit)[i], (*m_eLOC1_fit)[i], (*m_ePHI_fit)[i],
189  (*m_eTHETA_fit)[i], (*m_eQOP_fit)[i], (*m_eT_fit)[i];
190 
191  // Resolutions
192  double resD0 = (*m_err_eLOC0_fit)[i];
193  double resZ0 = (*m_err_eLOC1_fit)[i];
194  double resPh = (*m_err_ePHI_fit)[i];
195  double resTh = (*m_err_eTHETA_fit)[i];
196  double resQp = (*m_err_eQOP_fit)[i];
197  double resT = (*m_err_eT_fit)[i];
198 
199  // Fill vector of track objects with simple covariance matrix
200  Acts::BoundSquareMatrix covMat;
201 
202  covMat << resD0 * resD0, 0., 0., 0., 0., 0., 0., resZ0 * resZ0, 0., 0.,
203  0., 0., 0., 0., resPh * resPh, 0., 0., 0., 0., 0., 0., resTh * resTh,
204  0., 0., 0., 0., 0., 0., resQp * resQp, 0., 0., 0., 0., 0., 0.,
205  resT * resT;
206 
207  // TODO we do not have a hypothesis at hand here. defaulting to pion
208  trackParameterCollection.push_back(Acts::BoundTrackParameters(
209  perigeeSurface, paramVec, std::move(covMat),
210  Acts::ParticleHypothesis::pion()));
211  }
212 
213  unsigned int nTruthParticles = m_t_vx->size();
214  for (unsigned int i = 0; i < nTruthParticles; i++) {
215  ActsFatras::Particle truthParticle;
216 
217  truthParticle.setPosition4((*m_t_vx)[i], (*m_t_vy)[i], (*m_t_vz)[i],
218  (*m_t_time)[i]);
219  truthParticle.setDirection((*m_t_px)[i], (*m_t_py)[i], (*m_t_pz)[i]);
220  truthParticle.setParticleId((*m_majorityParticleId)[i]);
221 
222  truthParticleCollection.insert(truthParticleCollection.end(),
223  truthParticle);
224  }
225  // Write the collections to the EventStore
226  m_outputTrackParameters(context, std::move(trackParameterCollection));
227  m_outputParticles(context, std::move(truthParticleCollection));
228  } else {
229  ACTS_WARNING("Could not read in event.");
230  }
231  // Return success flag
232  return ActsExamples::ProcessCode::SUCCESS;
233 }