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