ECCE @ EIC Software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
RootParticleWriter.cpp
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file RootParticleWriter.cpp
1 // This file is part of the Acts project.
2 //
3 // Copyright (C) 2017 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 
10 
11 #include <TFile.h>
12 #include <TTree.h>
13 #include <ios>
14 #include <stdexcept>
15 
17 #include "Acts/Utilities/Units.hpp"
18 
21  : WriterT(cfg.inputParticles, "RootParticleWriter", 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("particle_id", &m_particleId, "particle_id/l");
44  m_outputTree->Branch("particle_type", &m_particleType);
45  m_outputTree->Branch("process", &m_process);
46  m_outputTree->Branch("vx", &m_vx);
47  m_outputTree->Branch("vy", &m_vy);
48  m_outputTree->Branch("vz", &m_vz);
49  m_outputTree->Branch("vt", &m_vt);
50  m_outputTree->Branch("px", &m_px);
51  m_outputTree->Branch("py", &m_py);
52  m_outputTree->Branch("pz", &m_pz);
53  m_outputTree->Branch("m", &m_m);
54  m_outputTree->Branch("q", &m_q);
55  m_outputTree->Branch("eta", &m_eta);
56  m_outputTree->Branch("phi", &m_phi);
57  m_outputTree->Branch("pt", &m_pt);
58  m_outputTree->Branch("vertex_primary", &m_vertexPrimary);
59  m_outputTree->Branch("vertex_secondary", &m_vertexSecondary);
60  m_outputTree->Branch("particle", &m_particle);
61  m_outputTree->Branch("generation", &m_generation);
62  m_outputTree->Branch("sub_particle", &m_subParticle);
63 }
64 
66  if (m_outputFile) {
67  m_outputFile->Close();
68  }
69 }
70 
72  if (m_outputFile) {
73  m_outputFile->cd();
74  m_outputTree->Write();
75  ACTS_INFO("Wrote particles to tree '" << m_cfg.treeName << "' in '"
76  << m_cfg.filePath << "'");
77  }
78  return ProcessCode::SUCCESS;
79 }
80 
82  const AlgorithmContext& ctx, const SimParticleContainer& particles) {
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  m_eventId = ctx.eventNumber;
92  for (const auto& particle : particles) {
93  m_particleId = particle.particleId().value();
94  m_particleType = particle.pdg();
95  m_process = static_cast<decltype(m_process)>(particle.process());
96  // position
97  m_vx = particle.position4().x() / Acts::UnitConstants::mm;
98  m_vy = particle.position4().y() / Acts::UnitConstants::mm;
99  m_vz = particle.position4().z() / Acts::UnitConstants::mm;
100  m_vt = particle.position4().w() / Acts::UnitConstants::ns;
101  // momentum
102  const auto p = particle.absMomentum() / Acts::UnitConstants::GeV;
103  m_px = p * particle.unitDirection().x();
104  m_py = p * particle.unitDirection().y();
105  m_pz = p * particle.unitDirection().z();
106  // particle constants
107  m_m = particle.mass() / Acts::UnitConstants::GeV;
108  m_q = particle.charge() / Acts::UnitConstants::e;
109  // derived kinematic quantities
110  m_eta = Acts::VectorHelpers::eta(particle.unitDirection());
111  m_phi = Acts::VectorHelpers::phi(particle.unitDirection());
112  m_pt = p * Acts::VectorHelpers::perp(particle.unitDirection());
113  // decoded barcode components
114  m_vertexPrimary = particle.particleId().vertexPrimary();
115  m_vertexSecondary = particle.particleId().vertexSecondary();
116  m_particle = particle.particleId().particle();
117  m_generation = particle.particleId().generation();
118  m_subParticle = particle.particleId().subParticle();
119  m_outputTree->Fill();
120  }
121 
122  return ProcessCode::SUCCESS;
123 }