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("barcode_vertex_primary", &m_barcodeVertexPrimary);
0051   m_outputTree->Branch("barcode_vertex_secondary", &m_barcodeVertexSecondary);
0052   m_outputTree->Branch("barcode_particle", &m_barcodeParticle);
0053   m_outputTree->Branch("barcode_generation", &m_barcodeGeneration);
0054   m_outputTree->Branch("barcode_sub_particle", &m_barcodeSubParticle);
0055   m_outputTree->Branch("tx", &m_tx);
0056   m_outputTree->Branch("ty", &m_ty);
0057   m_outputTree->Branch("tz", &m_tz);
0058   m_outputTree->Branch("tt", &m_tt);
0059   m_outputTree->Branch("tpx", &m_tpx);
0060   m_outputTree->Branch("tpy", &m_tpy);
0061   m_outputTree->Branch("tpz", &m_tpz);
0062   m_outputTree->Branch("te", &m_te);
0063   m_outputTree->Branch("deltapx", &m_deltapx);
0064   m_outputTree->Branch("deltapy", &m_deltapy);
0065   m_outputTree->Branch("deltapz", &m_deltapz);
0066   m_outputTree->Branch("deltae", &m_deltae);
0067   m_outputTree->Branch("index", &m_index);
0068   m_outputTree->Branch("volume_id", &m_volumeId);
0069   m_outputTree->Branch("boundary_id", &m_boundaryId);
0070   m_outputTree->Branch("layer_id", &m_layerId);
0071   m_outputTree->Branch("approach_id", &m_approachId);
0072   m_outputTree->Branch("sensitive_id", &m_sensitiveId);
0073 }
0074
0075 ActsExamples::RootSimHitWriter::~RootSimHitWriter() {
0076   if (m_outputFile != nullptr) {
0077     m_outputFile->Close();
0078   }
0079 }
0080
0081 ActsExamples::ProcessCode ActsExamples::RootSimHitWriter::finalize() {
0082   m_outputFile->cd();
0083   m_outputTree->Write();
0084   m_outputFile->Close();
0085
0086   ACTS_VERBOSE("Wrote hits to tree '" << m_cfg.treeName << "' in '"
0087                                       << m_cfg.filePath << "'");
0088
0089   return ProcessCode::SUCCESS;
0090 }
0091
0092 ActsExamples::ProcessCode ActsExamples::RootSimHitWriter::writeT(
0093     const AlgorithmContext& ctx, const ActsExamples::SimHitContainer& hits) {
0094   // ensure exclusive access to tree/file while writing
0095   std::lock_guard<std::mutex> lock(m_writeMutex);
0096
0097   // Get the event number
0098   m_eventId = ctx.eventNumber;
0099   for (const auto& hit : hits) {
0100     m_geometryId = hit.geometryId().value();
0101     const auto barcode = hit.particleId();
0102     m_barcodeVertexPrimary = barcode.vertexPrimary();
0103     m_barcodeVertexSecondary = barcode.vertexSecondary();
0104     m_barcodeParticle = barcode.particle();
0105     m_barcodeGeneration = barcode.generation();
0106     m_barcodeSubParticle = barcode.subParticle();
0107     // write hit position
0108     m_tx = hit.fourPosition().x() / Acts::UnitConstants::mm;
0109     m_ty = hit.fourPosition().y() / Acts::UnitConstants::mm;
0110     m_tz = hit.fourPosition().z() / Acts::UnitConstants::mm;
0111     m_tt = hit.fourPosition().w() / Acts::UnitConstants::mm;
0112     // write four-momentum before interaction
0113     m_tpx = hit.momentum4Before().x() / Acts::UnitConstants::GeV;
0114     m_tpy = hit.momentum4Before().y() / Acts::UnitConstants::GeV;
0115     m_tpz = hit.momentum4Before().z() / Acts::UnitConstants::GeV;
0116     m_te = hit.momentum4Before().w() / Acts::UnitConstants::GeV;
0117     // write four-momentum change due to interaction
0118     const auto delta4 = hit.momentum4After() - hit.momentum4Before();
0119     m_deltapx = delta4.x() / Acts::UnitConstants::GeV;
0120     m_deltapy = delta4.y() / Acts::UnitConstants::GeV;
0121     m_deltapz = delta4.z() / Acts::UnitConstants::GeV;
0122     m_deltae = delta4.w() / Acts::UnitConstants::GeV;
0123     // write hit index along trajectory
0124     m_index = hit.index();
0125     // decoded geometry for simplicity
0126     m_volumeId = hit.geometryId().volume();
0127     m_boundaryId = hit.geometryId().boundary();
0128     m_layerId = hit.geometryId().layer();
0129     m_approachId = hit.geometryId().approach();
0130     m_sensitiveId = hit.geometryId().sensitive();
0131     // Fill the tree
0132     m_outputTree->Fill();
0133   }
0134   return ActsExamples::ProcessCode::SUCCESS;
0135 }