Root/src/RootSimHitWriter.cpp

0001 // This file is part of the ACTS project.
0002 //
0003 // Copyright (C) 2016 CERN for the benefit of the ACTS project
0004 //
0005 // This Source Code Form is subject to the terms of the Mozilla Public
0006 // License, v. 2.0. If a copy of the MPL was not distributed with this
0007 // file, You can obtain one at https://mozilla.org/MPL/2.0/.
0008
0009 #include "ActsExamples/Io/Root/RootSimHitWriter.hpp"
0010
0011 #include "Acts/Definitions/Units.hpp"
0012 #include "Acts/Geometry/GeometryIdentifier.hpp"
0013 #include "ActsExamples/Framework/AlgorithmContext.hpp"
0014 #include "ActsFatras/EventData/Barcode.hpp"
0015 #include "ActsFatras/EventData/Hit.hpp"
0016
0017 #include <ios>
0018 #include <ostream>
0019 #include <stdexcept>
0020
0021 #include <TFile.h>
0022 #include <TTree.h>
0023
0024 ActsExamples::RootSimHitWriter::RootSimHitWriter(
0025     const ActsExamples::RootSimHitWriter::Config& config,
0026     Acts::Logging::Level level)
0027     : WriterT(config.inputSimHits, "RootSimHitWriter", level), m_cfg(config) {
0028   // inputParticles is already checked by base constructor
0029   if (m_cfg.filePath.empty()) {
0030     throw std::invalid_argument("Missing file path");
0031   }
0032   if (m_cfg.treeName.empty()) {
0033     throw std::invalid_argument("Missing tree name");
0034   }
0035
0036   // open root file and create the tree
0037   m_outputFile = TFile::Open(m_cfg.filePath.c_str(), m_cfg.fileMode.c_str());
0038   if (m_outputFile == nullptr) {
0039     throw std::ios_base::failure("Could not open '" + m_cfg.filePath + "'");
0040   }
0041   m_outputFile->cd();
0042   m_outputTree = new TTree(m_cfg.treeName.c_str(), m_cfg.treeName.c_str());
0043   if (m_outputTree == nullptr) {
0044     throw std::bad_alloc();
0045   }
0046
0047   // setup the branches
0048   m_outputTree->Branch("event_id", &m_eventId);
0049   m_outputTree->Branch("geometry_id", &m_geometryId, "geometry_id/l");
0050   m_outputTree->Branch("particle_id", &m_particleId, "particle_id/l");
0051   m_outputTree->Branch("tx", &m_tx);
0052   m_outputTree->Branch("ty", &m_ty);
0053   m_outputTree->Branch("tz", &m_tz);
0054   m_outputTree->Branch("tt", &m_tt);
0055   m_outputTree->Branch("tpx", &m_tpx);
0056   m_outputTree->Branch("tpy", &m_tpy);
0057   m_outputTree->Branch("tpz", &m_tpz);
0058   m_outputTree->Branch("te", &m_te);
0059   m_outputTree->Branch("deltapx", &m_deltapx);
0060   m_outputTree->Branch("deltapy", &m_deltapy);
0061   m_outputTree->Branch("deltapz", &m_deltapz);
0062   m_outputTree->Branch("deltae", &m_deltae);
0063   m_outputTree->Branch("index", &m_index);
0064   m_outputTree->Branch("volume_id", &m_volumeId);
0065   m_outputTree->Branch("boundary_id", &m_boundaryId);
0066   m_outputTree->Branch("layer_id", &m_layerId);
0067   m_outputTree->Branch("approach_id", &m_approachId);
0068   m_outputTree->Branch("sensitive_id", &m_sensitiveId);
0069 }
0070
0071 ActsExamples::RootSimHitWriter::~RootSimHitWriter() {
0072   if (m_outputFile != nullptr) {
0073     m_outputFile->Close();
0074   }
0075 }
0076
0077 ActsExamples::ProcessCode ActsExamples::RootSimHitWriter::finalize() {
0078   m_outputFile->cd();
0079   m_outputTree->Write();
0080   m_outputFile->Close();
0081
0082   ACTS_VERBOSE("Wrote hits to tree '" << m_cfg.treeName << "' in '"
0083                                       << m_cfg.filePath << "'");
0084
0085   return ProcessCode::SUCCESS;
0086 }
0087
0088 ActsExamples::ProcessCode ActsExamples::RootSimHitWriter::writeT(
0089     const AlgorithmContext& ctx, const ActsExamples::SimHitContainer& hits) {
0090   // ensure exclusive access to tree/file while writing
0091   std::lock_guard<std::mutex> lock(m_writeMutex);
0092
0093   // Get the event number
0094   m_eventId = ctx.eventNumber;
0095   for (const auto& hit : hits) {
0096     m_particleId = hit.particleId().value();
0097     m_geometryId = hit.geometryId().value();
0098     // write hit position
0099     m_tx = hit.fourPosition().x() / Acts::UnitConstants::mm;
0100     m_ty = hit.fourPosition().y() / Acts::UnitConstants::mm;
0101     m_tz = hit.fourPosition().z() / Acts::UnitConstants::mm;
0102     m_tt = hit.fourPosition().w() / Acts::UnitConstants::mm;
0103     // write four-momentum before interaction
0104     m_tpx = hit.momentum4Before().x() / Acts::UnitConstants::GeV;
0105     m_tpy = hit.momentum4Before().y() / Acts::UnitConstants::GeV;
0106     m_tpz = hit.momentum4Before().z() / Acts::UnitConstants::GeV;
0107     m_te = hit.momentum4Before().w() / Acts::UnitConstants::GeV;
0108     // write four-momentum change due to interaction
0109     const auto delta4 = hit.momentum4After() - hit.momentum4Before();
0110     m_deltapx = delta4.x() / Acts::UnitConstants::GeV;
0111     m_deltapy = delta4.y() / Acts::UnitConstants::GeV;
0112     m_deltapz = delta4.z() / Acts::UnitConstants::GeV;
0113     m_deltae = delta4.w() / Acts::UnitConstants::GeV;
0114     // write hit index along trajectory
0115     m_index = hit.index();
0116     // decoded geometry for simplicity
0117     m_volumeId = hit.geometryId().volume();
0118     m_boundaryId = hit.geometryId().boundary();
0119     m_layerId = hit.geometryId().layer();
0120     m_approachId = hit.geometryId().approach();
0121     m_sensitiveId = hit.geometryId().sensitive();
0122     // Fill the tree
0123     m_outputTree->Fill();
0124   }
0125   return ActsExamples::ProcessCode::SUCCESS;
0126 }