Root/src/RootPropagationSummaryWriter.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/RootPropagationSummaryWriter.hpp"
0010
0011 #include "Acts/Utilities/VectorHelpers.hpp"
0012 #include "ActsExamples/Framework/AlgorithmContext.hpp"
0013 #include "ActsExamples/Propagation/PropagationAlgorithm.hpp"
0014
0015 #include <ios>
0016 #include <ostream>
0017 #include <stdexcept>
0018
0019 #include <TFile.h>
0020 #include <TTree.h>
0021
0022 namespace ActsExamples {
0023
0024 RootPropagationSummaryWriter::RootPropagationSummaryWriter(
0025     const RootPropagationSummaryWriter::Config& cfg, Acts::Logging::Level level)
0026     : WriterT(cfg.inputSummaryCollection, "RootPropagationSummaryWriter",
0027               level),
0028       m_cfg(cfg),
0029       m_outputFile(cfg.rootFile) {
0030   // An input collection name and tree name must be specified
0031   if (m_cfg.inputSummaryCollection.empty()) {
0032     throw std::invalid_argument("Missing input collection");
0033   }
0034   if (m_cfg.treeName.empty()) {
0035     throw std::invalid_argument("Missing tree name");
0036   }
0037
0038   // Setup ROOT I/O
0039   if (m_outputFile == nullptr) {
0040     m_outputFile = TFile::Open(m_cfg.filePath.c_str(), m_cfg.fileMode.c_str());
0041     if (m_outputFile == nullptr) {
0042       throw std::ios_base::failure("Could not open '" + m_cfg.filePath + "'");
0043     }
0044   }
0045   m_outputFile->cd();
0046
0047   m_outputTree = new TTree(m_cfg.treeName.c_str(),
0048                            "TTree from RootPropagationSummaryWriter");
0049   if (m_outputTree == nullptr) {
0050     throw std::bad_alloc();
0051   }
0052
0053   // Set the branches
0054   m_outputTree->Branch("event_nr", &m_eventNr);
0055   m_outputTree->Branch("track_nr", &m_trackNr);
0056
0057   m_outputTree->Branch("d0", &m_d0);
0058   m_outputTree->Branch("z0", &m_z0);
0059   m_outputTree->Branch("phi", &m_phi);
0060   m_outputTree->Branch("theta", &m_theta);
0061   m_outputTree->Branch("qOverP", &m_qOverP);
0062   m_outputTree->Branch("t", &m_t);
0063
0064   m_outputTree->Branch("eta", &m_eta);
0065   m_outputTree->Branch("pt", &m_pt);
0066   m_outputTree->Branch("p", &m_p);
0067
0068   m_outputTree->Branch("nSensitives", &m_nSensitives);
0069   m_outputTree->Branch("nMaterials", &m_nMaterials);
0070   m_outputTree->Branch("nPortals", &m_nPortals);
0071
0072   m_outputTree->Branch("nAttemptedSteps", &m_nAttemptedSteps);
0073   m_outputTree->Branch("nRejectedSteps", &m_nRejectedSteps);
0074   m_outputTree->Branch("nSuccessfulSteps", &m_nSuccessfulSteps);
0075   m_outputTree->Branch("nReverseSteps", &m_nReverseSteps);
0076   m_outputTree->Branch("pathLength", &m_pathLength);
0077   m_outputTree->Branch("absolutePathLength", &m_absolutePathLength);
0078
0079   m_outputTree->Branch("nRenavigations", &m_nRenavigations);
0080   m_outputTree->Branch("nVolumeSwitches", &m_nVolumeSwitches);
0081 }
0082
0083 RootPropagationSummaryWriter::~RootPropagationSummaryWriter() {
0084   if (m_outputFile != nullptr) {
0085     m_outputFile->Close();
0086   }
0087 }
0088
0089 ProcessCode RootPropagationSummaryWriter::finalize() {
0090   // Write the tree
0091   m_outputFile->cd();
0092   m_outputTree->Write();
0093   /// Close the file if it's yours
0094   if (m_cfg.rootFile == nullptr) {
0095     m_outputFile->Close();
0096   }
0097
0098   ACTS_VERBOSE("Wrote particles to tree '" << m_cfg.treeName << "' in '"
0099                                            << m_cfg.filePath << "'");
0100
0101   return ProcessCode::SUCCESS;
0102 }
0103
0104 ProcessCode RootPropagationSummaryWriter::writeT(
0105     const AlgorithmContext& context, const PropagationSummaries& summaries) {
0106   // Exclusive access to the tree while writing
0107   std::scoped_lock lock(m_writeMutex);
0108
0109   // Get the event number
0110   m_eventNr = static_cast<int>(context.eventNumber);
0111
0112   // Loop over the step vector of each test propagation in this
0113   for (const auto& [trackNr, summary] : Acts::enumerate(summaries)) {
0114     // Set the track number
0115     m_trackNr = static_cast<int>(trackNr);
0116
0117     // Set the initial trajectory parameters
0118     const auto& startParameters = summary.startParameters;
0119     m_d0 = static_cast<float>(startParameters.parameters()[Acts::eBoundLoc0]);
0120     m_z0 = static_cast<float>(startParameters.parameters()[Acts::eBoundLoc1]);
0121     m_phi = static_cast<float>(startParameters.parameters()[Acts::eBoundPhi]);
0122     m_theta =
0123         static_cast<float>(startParameters.parameters()[Acts::eBoundTheta]);
0124     m_qOverP =
0125         static_cast<float>(startParameters.parameters()[Acts::eBoundQOverP]);
0126     m_t = static_cast<float>(startParameters.parameters()[Acts::eBoundTime]);
0127
0128     // Derived initial trajectory parameters
0129     m_eta = static_cast<float>(
0130         Acts::VectorHelpers::eta(startParameters.direction()));
0131     m_pt = static_cast<float>(startParameters.transverseMomentum());
0132     m_p = static_cast<float>(startParameters.absoluteMomentum());
0133
0134     m_nMaterials = 0;
0135     m_nSensitives = 0;
0136     m_nPortals = 0;
0137
0138     // Loop over the steps & count for the statistics
0139     std::ranges::for_each(summary.steps, [&](const auto& step) {
0140       // Check if the step is a sensitive step
0141       if (step.geoID.sensitive() > 0) {
0142         m_nSensitives++;
0143       }
0144       // Check if the step is a portal step
0145       if (step.geoID.boundary() > 0) {
0146         m_nPortals++;
0147       }
0148
0149       if (step.surface != nullptr) {
0150         // Check if the step is a material step
0151         if (step.surface->surfaceMaterial() != nullptr) {
0152           m_nMaterials++;
0153         }
0154       }
0155     });
0156
0157     // Stepper statistics
0158     m_nAttemptedSteps = summary.statistics.stepping.nAttemptedSteps;
0159     m_nRejectedSteps = summary.statistics.stepping.nRejectedSteps;
0160     m_nSuccessfulSteps = summary.statistics.stepping.nSuccessfulSteps;
0161     m_nReverseSteps = summary.statistics.stepping.nReverseSteps;
0162     m_pathLength = summary.statistics.stepping.pathLength;
0163     m_absolutePathLength = summary.statistics.stepping.absolutePathLength;
0164
0165     // Navigator statistics
0166     m_nRenavigations = summary.statistics.navigation.nRenavigations;
0167     m_nVolumeSwitches = summary.statistics.navigation.nVolumeSwitches;
0168
0169     m_outputTree->Fill();
0170   }
0171
0172   return ProcessCode::SUCCESS;
0173 }
0174
0175 }  // namespace ActsExamples