Analysis Software
Documentation for sPHENIX simulation software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
RootMaterialTrackReader.cpp
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file RootMaterialTrackReader.cpp
1 // This file is part of the Acts project.
2 //
3 // Copyright (C) 2017-2018 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/RootMaterialTrackReader.hpp"
10 
11 #include "Acts/Definitions/Algebra.hpp"
12 #include "Acts/Geometry/GeometryIdentifier.hpp"
13 #include "Acts/Material/Material.hpp"
14 #include "Acts/Material/MaterialSlab.hpp"
15 #include "Acts/Utilities/Logger.hpp"
16 #include "ActsExamples/Framework/AlgorithmContext.hpp"
17 
18 #include <cstdint>
19 #include <iostream>
20 #include <stdexcept>
21 
22 #include <TChain.h>
23 #include <TMathBase.h>
24 #include <TTree.h>
25 
26 ActsExamples::RootMaterialTrackReader::RootMaterialTrackReader(
27  const Config& config, Acts::Logging::Level level)
28  : ActsExamples::IReader(),
29  m_logger{Acts::getDefaultLogger(name(), level)},
30  m_cfg(config) {
31  if (m_cfg.fileList.empty()) {
32  throw std::invalid_argument{"No input files given"};
33  }
34 
35  m_inputChain = new TChain(m_cfg.treeName.c_str());
36 
37  // loop over the input files
38  for (const auto& inputFile : m_cfg.fileList) {
39  // add file to the input chain
40  m_inputChain->Add(inputFile.c_str());
41  ACTS_DEBUG("Adding File " << inputFile << " to tree '" << m_cfg.treeName
42  << "'.");
43  }
44 
45  // get the number of entries, which also loads the tree
46  size_t nentries = m_inputChain->GetEntries();
47 
48  bool eventIdPresent =
49  (TTree::kMatch == m_inputChain->SetBranchAddress("event_id", &m_eventId));
50 
51  m_inputChain->SetBranchAddress("v_x", &m_v_x);
52  m_inputChain->SetBranchAddress("v_y", &m_v_y);
53  m_inputChain->SetBranchAddress("v_z", &m_v_z);
54  m_inputChain->SetBranchAddress("v_px", &m_v_px);
55  m_inputChain->SetBranchAddress("v_py", &m_v_py);
56  m_inputChain->SetBranchAddress("v_pz", &m_v_pz);
57  m_inputChain->SetBranchAddress("v_phi", &m_v_phi);
58  m_inputChain->SetBranchAddress("v_eta", &m_v_eta);
59  m_inputChain->SetBranchAddress("t_X0", &m_tX0);
60  m_inputChain->SetBranchAddress("t_L0", &m_tL0);
61  m_inputChain->SetBranchAddress("mat_x", &m_step_x);
62  m_inputChain->SetBranchAddress("mat_y", &m_step_y);
63  m_inputChain->SetBranchAddress("mat_z", &m_step_z);
64  m_inputChain->SetBranchAddress("mat_dx", &m_step_dx);
65  m_inputChain->SetBranchAddress("mat_dy", &m_step_dy);
66  m_inputChain->SetBranchAddress("mat_dz", &m_step_dz);
67  m_inputChain->SetBranchAddress("mat_step_length", &m_step_length);
68  m_inputChain->SetBranchAddress("mat_X0", &m_step_X0);
69  m_inputChain->SetBranchAddress("mat_L0", &m_step_L0);
70  m_inputChain->SetBranchAddress("mat_A", &m_step_A);
71  m_inputChain->SetBranchAddress("mat_Z", &m_step_Z);
72  m_inputChain->SetBranchAddress("mat_rho", &m_step_rho);
73  if (m_cfg.readCachedSurfaceInformation) {
74  m_inputChain->SetBranchAddress("sur_id", &m_sur_id);
75  m_inputChain->SetBranchAddress("sur_x", &m_sur_x);
76  m_inputChain->SetBranchAddress("sur_y", &m_sur_y);
77  m_inputChain->SetBranchAddress("sur_z", &m_sur_z);
78  m_inputChain->SetBranchAddress("sur_pathCorrection", &m_sur_pathCorrection);
79  }
80 
81  m_events = eventIdPresent
82  ? static_cast<size_t>(m_inputChain->GetMaximum("event_id") + 1)
83  : nentries;
84  m_batchSize = nentries / m_events;
85  ACTS_DEBUG("The full chain has "
86  << nentries << " entries for " << m_events
87  << " events this corresponds to a batch size of: " << m_batchSize);
88  std::cout << "The full chain has " << nentries << " entries for " << m_events
89  << " events this corresponds to a batch size of: " << m_batchSize
90  << std::endl;
91 
92  // If the events are not in order, get the entry numbers for ordered events
93  if (not m_cfg.orderedEvents) {
94  if (not eventIdPresent) {
95  throw std::invalid_argument{
96  "'event_id' branch is missing in your tree. This is not compatible "
97  "with unordered events."};
98  }
99  m_entryNumbers.resize(nentries);
100  m_inputChain->Draw("event_id", "", "goff");
101  // Sort to get the entry numbers of the ordered events
102  TMath::Sort(m_inputChain->GetEntries(), m_inputChain->GetV1(),
103  m_entryNumbers.data(), false);
104  }
105  m_outputMaterialTracks.initialize(m_cfg.collection);
106 }
107 
108 ActsExamples::RootMaterialTrackReader::~RootMaterialTrackReader() {
109  delete m_inputChain;
110 
111  delete m_step_x;
112  delete m_step_y;
113  delete m_step_z;
114  delete m_step_dx;
115  delete m_step_dy;
116  delete m_step_dz;
117  delete m_step_length;
118  delete m_step_X0;
119  delete m_step_L0;
120  delete m_step_A;
121  delete m_step_Z;
122  delete m_step_rho;
123 
124  delete m_sur_id;
125  delete m_sur_x;
126  delete m_sur_y;
127  delete m_sur_z;
128  delete m_sur_pathCorrection;
129 }
130 
131 std::string ActsExamples::RootMaterialTrackReader::name() const {
132  return "RootMaterialTrackReader";
133 }
134 
135 std::pair<size_t, size_t>
136 ActsExamples::RootMaterialTrackReader::availableEvents() const {
137  return {0u, m_events};
138 }
139 
140 ActsExamples::ProcessCode ActsExamples::RootMaterialTrackReader::read(
141  const ActsExamples::AlgorithmContext& context) {
142  ACTS_DEBUG("Trying to read recorded material from tracks.");
143  // read in the material track
144  if (m_inputChain != nullptr && context.eventNumber < m_events) {
145  // lock the mutex
146  std::lock_guard<std::mutex> lock(m_read_mutex);
147  // now read
148 
149  // The collection to be written
150  std::unordered_map<size_t, Acts::RecordedMaterialTrack> mtrackCollection;
151 
152  // Loop over the entries for this event
153  for (size_t ib = 0; ib < m_batchSize; ++ib) {
154  // Read the correct entry: startEntry + ib
155  auto entry = m_batchSize * context.eventNumber + ib;
156  if (not m_cfg.orderedEvents and entry < m_entryNumbers.size()) {
157  entry = m_entryNumbers[entry];
158  }
159  ACTS_VERBOSE("Reading event: " << context.eventNumber
160  << " with stored entry: " << entry);
161  m_inputChain->GetEntry(entry);
162 
163  Acts::RecordedMaterialTrack rmTrack;
164  // Fill the position and momentum
165  rmTrack.first.first = Acts::Vector3(m_v_x, m_v_y, m_v_z);
166  rmTrack.first.second = Acts::Vector3(m_v_px, m_v_py, m_v_pz);
167 
168  ACTS_VERBOSE("Track vertex: " << rmTrack.first.first);
169  ACTS_VERBOSE("Track momentum:" << rmTrack.first.second);
170 
171  // Fill the individual steps
172  size_t msteps = m_step_length->size();
173  ACTS_VERBOSE("Reading " << msteps << " material steps.");
174  rmTrack.second.materialInteractions.reserve(msteps);
175  rmTrack.second.materialInX0 = 0.;
176  rmTrack.second.materialInL0 = 0.;
177 
178  for (size_t is = 0; is < msteps; ++is) {
179  ACTS_VERBOSE("====================");
180  ACTS_VERBOSE("[" << is + 1 << "/" << msteps << "] STEP INFORMATION: ");
181 
182  double s = (*m_step_length)[is];
183  if (s == 0) {
184  ACTS_VERBOSE("invalid step length... skipping!");
185  continue;
186  }
187 
188  double mX0 = (*m_step_X0)[is];
189  double mL0 = (*m_step_L0)[is];
190 
191  rmTrack.second.materialInX0 += s / mX0;
192  rmTrack.second.materialInL0 += s / mL0;
194  Acts::MaterialInteraction mInteraction;
195  mInteraction.position =
196  Acts::Vector3((*m_step_x)[is], (*m_step_y)[is], (*m_step_z)[is]);
197  ACTS_VERBOSE("POSITION : " << (*m_step_x)[is] << ", " << (*m_step_y)[is]
198  << ", " << (*m_step_z)[is]);
199  mInteraction.direction =
200  Acts::Vector3((*m_step_dx)[is], (*m_step_dy)[is], (*m_step_dz)[is]);
201  ACTS_VERBOSE("DIRECTION: " << (*m_step_dx)[is] << ", "
202  << (*m_step_dy)[is] << ", "
203  << (*m_step_dz)[is]);
204  mInteraction.materialSlab = Acts::MaterialSlab(
205  Acts::Material::fromMassDensity(mX0, mL0, (*m_step_A)[is],
206  (*m_step_Z)[is], (*m_step_rho)[is]),
207  s);
208  ACTS_VERBOSE("MATERIAL: " << mX0 << ", " << mL0 << ", "
209  << (*m_step_A)[is] << ", " << (*m_step_Z)[is]
210  << ", " << (*m_step_rho)[is]);
211  ACTS_VERBOSE("====================");
212 
213  if (m_cfg.readCachedSurfaceInformation) {
214  // add the surface information to the interaction this allows the
215  // mapping to be speed up
216  mInteraction.intersectionID =
217  Acts::GeometryIdentifier((*m_sur_id)[is]);
218  mInteraction.intersection =
219  Acts::Vector3((*m_sur_x)[is], (*m_sur_y)[is], (*m_sur_z)[is]);
220  mInteraction.pathCorrection = (*m_sur_pathCorrection)[is];
221  } else {
222  mInteraction.intersectionID = Acts::GeometryIdentifier();
223  mInteraction.intersection = Acts::Vector3(0, 0, 0);
224  }
225  rmTrack.second.materialInteractions.push_back(std::move(mInteraction));
226  }
227  mtrackCollection[ib] = (std::move(rmTrack));
228  }
229  // Write to the collection to the EventStore
230  m_outputMaterialTracks(context, std::move(mtrackCollection));
231  }
232  // Return success flag
233  return ActsExamples::ProcessCode::SUCCESS;
234 }