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