EIC Software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
RootSimHitWriter.cpp
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file RootSimHitWriter.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/RootSimHitWriter.hpp"
10 
11 #include "Acts/Utilities/Units.hpp"
12 
13 #include <ios>
14 #include <stdexcept>
15 
16 #include <TFile.h>
17 #include <TTree.h>
18 
19 ActsExamples::RootSimHitWriter::RootSimHitWriter(
20  const ActsExamples::RootSimHitWriter::Config& cfg, Acts::Logging::Level lvl)
21  : WriterT(cfg.inputSimulatedHits, "RootSimHitWriter", lvl), m_cfg(cfg) {
22  // inputParticles is already checked by base constructor
23  if (m_cfg.filePath.empty()) {
24  throw std::invalid_argument("Missing file path");
25  }
26  if (m_cfg.treeName.empty()) {
27  throw std::invalid_argument("Missing tree name");
28  }
29 
30  // open root file and create the tree
31  m_outputFile = TFile::Open(m_cfg.filePath.c_str(), m_cfg.fileMode.c_str());
32  if (m_outputFile == nullptr) {
33  throw std::ios_base::failure("Could not open '" + m_cfg.filePath + "'");
34  }
35  m_outputFile->cd();
36  m_outputTree = new TTree(m_cfg.treeName.c_str(), m_cfg.treeName.c_str());
37  if (m_outputTree == nullptr) {
38  throw std::bad_alloc();
39  }
40 
41  // setup the branches
42  m_outputTree->Branch("event_id", &m_eventId);
43  m_outputTree->Branch("geometry_id", &m_geometryId, "geometry_id/l");
44  m_outputTree->Branch("particle_id", &m_particleId, "particle_id/l");
45  m_outputTree->Branch("tx", &m_tx);
46  m_outputTree->Branch("ty", &m_ty);
47  m_outputTree->Branch("tz", &m_tz);
48  m_outputTree->Branch("tt", &m_tt);
49  m_outputTree->Branch("tpx", &m_tpx);
50  m_outputTree->Branch("tpy", &m_tpy);
51  m_outputTree->Branch("tpz", &m_tpz);
52  m_outputTree->Branch("te", &m_te);
53  m_outputTree->Branch("deltapx", &m_deltapx);
54  m_outputTree->Branch("deltapy", &m_deltapy);
55  m_outputTree->Branch("deltapz", &m_deltapz);
56  m_outputTree->Branch("deltae", &m_deltae);
57  m_outputTree->Branch("index", &m_index);
58  m_outputTree->Branch("volume_id", &m_volumeId);
59  m_outputTree->Branch("boundary_id", &m_boundaryId);
60  m_outputTree->Branch("layer_id", &m_layerId);
61  m_outputTree->Branch("approach_id", &m_approachId);
62  m_outputTree->Branch("sensitive_id", &m_sensitiveId);
63 }
64 
65 ActsExamples::RootSimHitWriter::~RootSimHitWriter() {
66  if (m_outputFile) {
67  m_outputFile->Close();
68  }
69 }
70 
71 ActsExamples::ProcessCode ActsExamples::RootSimHitWriter::endRun() {
72  if (m_outputFile) {
73  m_outputFile->cd();
74  m_outputTree->Write();
75  ACTS_VERBOSE("Wrote hits to tree '" << m_cfg.treeName << "' in '"
76  << m_cfg.filePath << "'");
77  }
78  return ProcessCode::SUCCESS;
79 }
80 
81 ActsExamples::ProcessCode ActsExamples::RootSimHitWriter::writeT(
82  const AlgorithmContext& ctx, const ActsExamples::SimHitContainer& hits) {
83  if (not m_outputFile) {
84  ACTS_ERROR("Missing output file");
85  return ProcessCode::ABORT;
86  }
87 
88  // ensure exclusive access to tree/file while writing
89  std::lock_guard<std::mutex> lock(m_writeMutex);
90 
91  // Get the event number
92  m_eventId = ctx.eventNumber;
93  for (const auto& hit : hits) {
94  m_particleId = hit.particleId().value();
95  m_geometryId = hit.geometryId().value();
96  // write hit position
97  m_tx = hit.position4().x() / Acts::UnitConstants::mm;
98  m_ty = hit.position4().y() / Acts::UnitConstants::mm;
99  m_tz = hit.position4().z() / Acts::UnitConstants::mm;
100  m_tt = hit.position4().w() / Acts::UnitConstants::ns;
101  // write four-momentum before interaction
102  m_tpx = hit.momentum4Before().x() / Acts::UnitConstants::GeV;
103  m_tpy = hit.momentum4Before().y() / Acts::UnitConstants::GeV;
104  m_tpz = hit.momentum4Before().z() / Acts::UnitConstants::GeV;
105  m_te = hit.momentum4Before().w() / Acts::UnitConstants::GeV;
106  // write four-momentum change due to interaction
107  const auto delta4 = hit.momentum4After() - hit.momentum4Before();
108  m_deltapx = delta4.x() / Acts::UnitConstants::GeV;
109  m_deltapy = delta4.y() / Acts::UnitConstants::GeV;
110  m_deltapz = delta4.z() / Acts::UnitConstants::GeV;
111  m_deltae = delta4.w() / Acts::UnitConstants::GeV;
112  // write hit index along trajectory
113  m_index = hit.index();
114  // decoded geometry for simplicity
115  m_volumeId = hit.geometryId().volume();
116  m_boundaryId = hit.geometryId().boundary();
117  m_layerId = hit.geometryId().layer();
118  m_approachId = hit.geometryId().approach();
119  m_sensitiveId = hit.geometryId().sensitive();
120  // Fill the tree
121  m_outputTree->Fill();
122  }
123  return ActsExamples::ProcessCode::SUCCESS;
124 }