ECCE @ EIC Software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
RootMaterialTrackWriter.cpp
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file RootMaterialTrackWriter.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 "ACTFW/Io/Root/RootMaterialTrackWriter.hpp"
10 
11 #include <Acts/Geometry/GeometryID.hpp>
12 #include <Acts/Surfaces/CylinderBounds.hpp>
13 #include <Acts/Surfaces/RadialBounds.hpp>
14 #include <Acts/Utilities/Helpers.hpp>
15 #include <TFile.h>
16 #include <TTree.h>
17 #include <ios>
18 #include <iostream>
19 #include <stdexcept>
20 
21 using Acts::VectorHelpers::eta;
22 using Acts::VectorHelpers::perp;
23 using Acts::VectorHelpers::phi;
24 
25 FW::RootMaterialTrackWriter::RootMaterialTrackWriter(
26  const FW::RootMaterialTrackWriter::Config& cfg, Acts::Logging::Level level)
27  : WriterT(cfg.collection, "RootMaterialTrackWriter", level),
28  m_cfg(cfg),
29  m_outputFile(cfg.rootFile) {
30  // An input collection name and tree name must be specified
31  if (m_cfg.collection.empty()) {
32  throw std::invalid_argument("Missing input collection");
33  } else if (m_cfg.treeName.empty()) {
34  throw std::invalid_argument("Missing tree name");
35  }
36 
37  // Setup ROOT I/O
38  if (m_outputFile == nullptr) {
39  m_outputFile = TFile::Open(m_cfg.filePath.c_str(), m_cfg.fileMode.c_str());
40  if (m_outputFile == nullptr) {
41  throw std::ios_base::failure("Could not open '" + m_cfg.filePath);
42  }
43  }
44  m_outputFile->cd();
45  m_outputTree =
46  new TTree(m_cfg.treeName.c_str(), "TTree from RootMaterialTrackWriter");
47  if (m_outputTree == nullptr)
48  throw std::bad_alloc();
49 
50  // Set the branches
51  m_outputTree->Branch("v_x", &m_v_x);
52  m_outputTree->Branch("v_y", &m_v_y);
53  m_outputTree->Branch("v_z", &m_v_z);
54  m_outputTree->Branch("v_px", &m_v_px);
55  m_outputTree->Branch("v_py", &m_v_py);
56  m_outputTree->Branch("v_pz", &m_v_pz);
57  m_outputTree->Branch("v_phi", &m_v_phi);
58  m_outputTree->Branch("v_eta", &m_v_eta);
59  m_outputTree->Branch("t_X0", &m_tX0);
60  m_outputTree->Branch("t_L0", &m_tL0);
61  m_outputTree->Branch("mat_x", &m_step_x);
62  m_outputTree->Branch("mat_y", &m_step_y);
63  m_outputTree->Branch("mat_z", &m_step_z);
64  m_outputTree->Branch("mat_step_length", &m_step_length);
65  m_outputTree->Branch("mat_X0", &m_step_X0);
66  m_outputTree->Branch("mat_L0", &m_step_L0);
67  m_outputTree->Branch("mat_A", &m_step_A);
68  m_outputTree->Branch("mat_Z", &m_step_Z);
69  m_outputTree->Branch("mat_rho", &m_step_rho);
70 
71  if (m_cfg.prePostStep) {
72  m_outputTree->Branch("mat_sx", &m_step_sx);
73  m_outputTree->Branch("mat_sy", &m_step_sy);
74  m_outputTree->Branch("mat_sz", &m_step_sz);
75  m_outputTree->Branch("mat_ex", &m_step_ex);
76  m_outputTree->Branch("mat_ey", &m_step_ey);
77  m_outputTree->Branch("mat_ez", &m_step_ez);
78  }
79  if (m_cfg.storesurface) {
80  m_outputTree->Branch("sur_id", &m_sur_id);
81  m_outputTree->Branch("sur_type", &m_sur_type);
82  m_outputTree->Branch("sur_x", &m_sur_x);
83  m_outputTree->Branch("sur_y", &m_sur_y);
84  m_outputTree->Branch("sur_z", &m_sur_z);
85  m_outputTree->Branch("sur_range_min", &m_sur_range_min);
86  m_outputTree->Branch("sur_range_max", &m_sur_range_max);
87  }
88 }
89 
90 FW::RootMaterialTrackWriter::~RootMaterialTrackWriter() {
91  m_outputFile->Close();
92 }
93 
94 FW::ProcessCode FW::RootMaterialTrackWriter::endRun() {
95  // write the tree and close the file
96  ACTS_INFO("Writing ROOT output File : " << m_cfg.filePath);
97  m_outputFile->cd();
98  m_outputTree->Write();
99  return FW::ProcessCode::SUCCESS;
100 }
101 
102 FW::ProcessCode FW::RootMaterialTrackWriter::writeT(
103  const AlgorithmContext& ctx,
104  const std::vector<Acts::RecordedMaterialTrack>& materialTracks) {
105  // Exclusive access to the tree while writing
106  std::lock_guard<std::mutex> lock(m_writeMutex);
107 
108  // Loop over the material tracks and write them out
109  for (auto& mtrack : materialTracks) {
110  // Clearing the vector first
111  m_step_sx.clear();
112  m_step_sy.clear();
113  m_step_sz.clear();
114  m_step_x.clear();
115  m_step_y.clear();
116  m_step_z.clear();
117  m_step_ex.clear();
118  m_step_ey.clear();
119  m_step_ez.clear();
120  m_step_length.clear();
121  m_step_X0.clear();
122  m_step_L0.clear();
123  m_step_A.clear();
124  m_step_Z.clear();
125  m_step_rho.clear();
126 
127  m_sur_id.clear();
128  m_sur_type.clear();
129  m_sur_x.clear();
130  m_sur_y.clear();
131  m_sur_z.clear();
132  m_sur_range_min.clear();
133  m_sur_range_max.clear();
134 
135  // Reserve the vector then
136  size_t mints = mtrack.second.materialInteractions.size();
137  m_step_sx.reserve(mints);
138  m_step_sy.reserve(mints);
139  m_step_sz.reserve(mints);
140  m_step_x.reserve(mints);
141  m_step_y.reserve(mints);
142  m_step_ez.reserve(mints);
143  m_step_ex.reserve(mints);
144  m_step_ey.reserve(mints);
145  m_step_ez.reserve(mints);
146  m_step_length.reserve(mints);
147  m_step_X0.reserve(mints);
148  m_step_L0.reserve(mints);
149  m_step_A.reserve(mints);
150  m_step_Z.reserve(mints);
151  m_step_rho.reserve(mints);
152 
153  m_sur_id.reserve(mints);
154  m_sur_type.reserve(mints);
155  m_sur_x.reserve(mints);
156  m_sur_y.reserve(mints);
157  m_sur_z.reserve(mints);
158  m_sur_range_min.reserve(mints);
159  m_sur_range_max.reserve(mints);
160 
161  // reset the global counter
162  if (m_cfg.recalculateTotals) {
163  m_tX0 = 0.;
164  m_tL0 = 0.;
165  } else {
166  m_tX0 = mtrack.second.materialInX0;
167  m_tL0 = mtrack.second.materialInL0;
168  }
169 
170  // set the track information at vertex
171  m_v_x = mtrack.first.first.x();
172  m_v_y = mtrack.first.first.y();
173  m_v_z = mtrack.first.first.z();
174  m_v_px = mtrack.first.second.x();
175  m_v_py = mtrack.first.second.y();
176  m_v_pz = mtrack.first.second.z();
177  m_v_phi = phi(mtrack.first.second);
178  m_v_eta = eta(mtrack.first.second);
179 
180  // an now loop over the material
181  for (auto& mint : mtrack.second.materialInteractions) {
182  // The material step position information
183  m_step_x.push_back(mint.position.x());
184  m_step_y.push_back(mint.position.y());
185  m_step_z.push_back(mint.position.z());
186 
187  if (m_cfg.prePostStep) {
188  Acts::Vector3D prePos =
189  mint.position - 0.5 * mint.pathCorrection * mint.direction;
190  Acts::Vector3D posPos =
191  mint.position + 0.5 * mint.pathCorrection * mint.direction;
192  m_step_sx.push_back(prePos.x());
193  m_step_sy.push_back(prePos.y());
194  m_step_sz.push_back(prePos.z());
195  m_step_ex.push_back(posPos.x());
196  m_step_ey.push_back(posPos.y());
197  m_step_ez.push_back(posPos.z());
198  }
199 
200  if (m_cfg.storesurface) {
201  const Acts::Surface* surface = mint.surface;
202  Acts::GeometryID layerID;
203  if (surface) {
204  Acts::Intersection intersection = surface->intersectionEstimate(
205  ctx.geoContext, mint.position, mint.direction, true);
206  layerID = surface->geoID();
207  m_sur_id.push_back(layerID.value());
208  m_sur_type.push_back(surface->type());
209  m_sur_x.push_back(intersection.position.x());
210  m_sur_y.push_back(intersection.position.y());
211  m_sur_z.push_back(intersection.position.z());
212 
213  const Acts::SurfaceBounds& surfaceBounds = surface->bounds();
214  const Acts::RadialBounds* radialBounds =
215  dynamic_cast<const Acts::RadialBounds*>(&surfaceBounds);
216  const Acts::CylinderBounds* cylinderBounds =
217  dynamic_cast<const Acts::CylinderBounds*>(&surfaceBounds);
218 
219  if (radialBounds) {
220  m_sur_range_min.push_back(radialBounds->rMin());
221  m_sur_range_max.push_back(radialBounds->rMax());
222  } else if (cylinderBounds) {
223  m_sur_range_min.push_back(
224  -cylinderBounds->get(Acts::CylinderBounds::eHalfLengthZ));
225  m_sur_range_max.push_back(
226  cylinderBounds->get(Acts::CylinderBounds::eHalfLengthZ));
227  } else {
228  m_sur_range_min.push_back(0);
229  m_sur_range_max.push_back(0);
230  }
231  } else {
232  layerID.setVolume(0);
233  layerID.setBoundary(0);
234  layerID.setLayer(0);
235  layerID.setApproach(0);
236  layerID.setSensitive(0);
237  m_sur_id.push_back(layerID.value());
238  m_sur_type.push_back(-1);
239 
240  m_sur_x.push_back(0);
241  m_sur_y.push_back(0);
242  m_sur_z.push_back(0);
243  m_sur_range_min.push_back(0);
244  m_sur_range_max.push_back(0);
245  }
246  }
247 
248  // the material information
249  const auto& mprops = mint.materialProperties;
250  m_step_length.push_back(mprops.thickness());
251  m_step_X0.push_back(mprops.material().X0());
252  m_step_L0.push_back(mprops.material().L0());
253  m_step_A.push_back(mprops.material().Ar());
254  m_step_Z.push_back(mprops.material().Z());
255  m_step_rho.push_back(mprops.material().massDensity());
256  // re-calculate if defined to do so
257  if (m_cfg.recalculateTotals) {
258  m_tX0 += mprops.thicknessInX0();
259  m_tL0 += mprops.thicknessInL0();
260  }
261  }
262  // write to
263  m_outputTree->Fill();
264  }
265 
266  // return success
267  return FW::ProcessCode::SUCCESS;
268 }