Root/src/RootMaterialTrackReader.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/RootMaterialTrackReader.hpp"
0010
0011 #include "Acts/Definitions/Algebra.hpp"
0012 #include "Acts/Geometry/GeometryIdentifier.hpp"
0013 #include "Acts/Material/Material.hpp"
0014 #include "Acts/Material/MaterialSlab.hpp"
0015 #include "Acts/Utilities/Logger.hpp"
0016 #include "ActsExamples/Framework/AlgorithmContext.hpp"
0017 #include "ActsExamples/Io/Root/RootUtility.hpp"
0018
0019 #include <cstdint>
0020 #include <iostream>
0021 #include <stdexcept>
0022
0023 #include <TChain.h>
0024 #include <TTree.h>
0025
0026 namespace ActsExamples {
0027
0028 RootMaterialTrackReader::RootMaterialTrackReader(const Config& config,
0029                                                  Acts::Logging::Level level)
0030     : IReader(),
0031       m_logger{Acts::getDefaultLogger(name(), level)},
0032       m_cfg(config) {
0033   if (m_cfg.fileList.empty()) {
0034     throw std::invalid_argument{"No input files given"};
0035   }
0036
0037   m_inputChain = std::make_unique<TChain>(m_cfg.treeName.c_str());
0038
0039   // loop over the input files
0040   for (const auto& inputFile : m_cfg.fileList) {
0041     // add file to the input chain
0042     m_inputChain->Add(inputFile.c_str());
0043     ACTS_DEBUG("Adding File " << inputFile << " to tree '" << m_cfg.treeName
0044                               << "'.");
0045   }
0046
0047   // get the number of entries, which also loads the tree
0048   std::size_t nentries = m_inputChain->GetEntries();
0049
0050   m_inputChain->SetBranchAddress("event_id", &m_eventId);
0051   m_inputChain->SetBranchAddress("v_x", &m_v_x);
0052   m_inputChain->SetBranchAddress("v_y", &m_v_y);
0053   m_inputChain->SetBranchAddress("v_z", &m_v_z);
0054   m_inputChain->SetBranchAddress("v_px", &m_v_px);
0055   m_inputChain->SetBranchAddress("v_py", &m_v_py);
0056   m_inputChain->SetBranchAddress("v_pz", &m_v_pz);
0057   m_inputChain->SetBranchAddress("v_phi", &m_v_phi);
0058   m_inputChain->SetBranchAddress("v_eta", &m_v_eta);
0059   m_inputChain->SetBranchAddress("t_X0", &m_tX0);
0060   m_inputChain->SetBranchAddress("t_L0", &m_tL0);
0061   m_inputChain->SetBranchAddress("mat_x", &m_step_x);
0062   m_inputChain->SetBranchAddress("mat_y", &m_step_y);
0063   m_inputChain->SetBranchAddress("mat_z", &m_step_z);
0064   m_inputChain->SetBranchAddress("mat_dx", &m_step_dx);
0065   m_inputChain->SetBranchAddress("mat_dy", &m_step_dy);
0066   m_inputChain->SetBranchAddress("mat_dz", &m_step_dz);
0067   m_inputChain->SetBranchAddress("mat_step_length", &m_step_length);
0068   m_inputChain->SetBranchAddress("mat_X0", &m_step_X0);
0069   m_inputChain->SetBranchAddress("mat_L0", &m_step_L0);
0070   m_inputChain->SetBranchAddress("mat_A", &m_step_A);
0071   m_inputChain->SetBranchAddress("mat_Z", &m_step_Z);
0072   m_inputChain->SetBranchAddress("mat_rho", &m_step_rho);
0073   if (m_cfg.readCachedSurfaceInformation) {
0074     m_inputChain->SetBranchAddress("sur_id", &m_sur_id);
0075     m_inputChain->SetBranchAddress("sur_x", &m_sur_x);
0076     m_inputChain->SetBranchAddress("sur_y", &m_sur_y);
0077     m_inputChain->SetBranchAddress("sur_z", &m_sur_z);
0078     m_inputChain->SetBranchAddress("sur_pathCorrection", &m_sur_pathCorrection);
0079   }
0080
0081   m_events = static_cast<std::size_t>(m_inputChain->GetMaximum("event_id") + 1);
0082   m_batchSize = nentries / m_events;
0083   ACTS_DEBUG("The full chain has "
0084              << nentries << " entries for " << m_events
0085              << " events this corresponds to a batch size of: " << m_batchSize);
0086
0087   // Sort the entry numbers of the events
0088   {
0089     // necessary to guarantee that m_inputChain->GetV1() is valid for the
0090     // entire range
0091     m_inputChain->SetEstimate(nentries + 1);
0092
0093     m_entryNumbers.resize(nentries);
0094     m_inputChain->Draw("event_id", "", "goff");
0095     RootUtility::stableSort(m_inputChain->GetEntries(), m_inputChain->GetV1(),
0096                             m_entryNumbers.data(), false);
0097   }
0098
0099   m_outputMaterialTracks.initialize(m_cfg.outputMaterialTracks);
0100 }
0101
0102 RootMaterialTrackReader::~RootMaterialTrackReader() {
0103   delete m_step_x;
0104   delete m_step_y;
0105   delete m_step_z;
0106   delete m_step_dx;
0107   delete m_step_dy;
0108   delete m_step_dz;
0109   delete m_step_length;
0110   delete m_step_X0;
0111   delete m_step_L0;
0112   delete m_step_A;
0113   delete m_step_Z;
0114   delete m_step_rho;
0115
0116   delete m_sur_id;
0117   delete m_sur_x;
0118   delete m_sur_y;
0119   delete m_sur_z;
0120   delete m_sur_pathCorrection;
0121 }
0122
0123 std::string RootMaterialTrackReader::name() const {
0124   return "RootMaterialTrackReader";
0125 }
0126
0127 std::pair<std::size_t, std::size_t> RootMaterialTrackReader::availableEvents()
0128     const {
0129   return {0u, m_events};
0130 }
0131
0132 ProcessCode RootMaterialTrackReader::read(const AlgorithmContext& context) {
0133   ACTS_DEBUG("Trying to read recorded material from tracks.");
0134
0135   if (m_inputChain == nullptr || context.eventNumber >= m_events) {
0136     return ProcessCode::SUCCESS;
0137   }
0138
0139   // lock the mutex
0140   std::lock_guard<std::mutex> lock(m_read_mutex);
0141   // now read
0142
0143   // The collection to be written
0144   std::unordered_map<std::size_t, Acts::RecordedMaterialTrack> mtrackCollection;
0145
0146   // Loop over the entries for this event
0147   for (std::size_t ib = 0; ib < m_batchSize; ++ib) {
0148     // Read the correct entry: startEntry + ib
0149     auto entry = m_batchSize * context.eventNumber + ib;
0150     entry = m_entryNumbers.at(entry);
0151     ACTS_VERBOSE("Reading event: " << context.eventNumber
0152                                    << " with stored entry: " << entry);
0153     m_inputChain->GetEntry(entry);
0154
0155     Acts::RecordedMaterialTrack rmTrack;
0156     // Fill the position and momentum
0157     rmTrack.first.first = Acts::Vector3(m_v_x, m_v_y, m_v_z);
0158     rmTrack.first.second = Acts::Vector3(m_v_px, m_v_py, m_v_pz);
0159
0160     ACTS_VERBOSE("Track vertex:  " << rmTrack.first.first);
0161     ACTS_VERBOSE("Track momentum:" << rmTrack.first.second);
0162
0163     // Fill the individual steps
0164     std::size_t msteps = m_step_length->size();
0165     ACTS_VERBOSE("Reading " << msteps << " material steps.");
0166     rmTrack.second.materialInteractions.reserve(msteps);
0167     rmTrack.second.materialInX0 = 0.;
0168     rmTrack.second.materialInL0 = 0.;
0169
0170     for (std::size_t is = 0; is < msteps; ++is) {
0171       ACTS_VERBOSE("====================");
0172       ACTS_VERBOSE("[" << is + 1 << "/" << msteps << "] STEP INFORMATION: ");
0173
0174       double s = (*m_step_length)[is];
0175       if (s == 0) {
0176         ACTS_VERBOSE("invalid step length... skipping!");
0177         continue;
0178       }
0179
0180       double mX0 = (*m_step_X0)[is];
0181       double mL0 = (*m_step_L0)[is];
0182
0183       rmTrack.second.materialInX0 += s / mX0;
0184       rmTrack.second.materialInL0 += s / mL0;
0185       /// Fill the position & the material
0186       Acts::MaterialInteraction mInteraction;
0187       mInteraction.position =
0188           Acts::Vector3((*m_step_x)[is], (*m_step_y)[is], (*m_step_z)[is]);
0189       ACTS_VERBOSE("POSITION : " << (*m_step_x)[is] << ", " << (*m_step_y)[is]
0190                                  << ", " << (*m_step_z)[is]);
0191       mInteraction.direction =
0192           Acts::Vector3((*m_step_dx)[is], (*m_step_dy)[is], (*m_step_dz)[is]);
0193       ACTS_VERBOSE("DIRECTION: " << (*m_step_dx)[is] << ", " << (*m_step_dy)[is]
0194                                  << ", " << (*m_step_dz)[is]);
0195       mInteraction.materialSlab = Acts::MaterialSlab(
0196           Acts::Material::fromMassDensity(mX0, mL0, (*m_step_A)[is],
0197                                           (*m_step_Z)[is], (*m_step_rho)[is]),
0198           s);
0199       ACTS_VERBOSE("MATERIAL: " << mX0 << ", " << mL0 << ", " << (*m_step_A)[is]
0200                                 << ", " << (*m_step_Z)[is] << ", "
0201                                 << (*m_step_rho)[is]);
0202       ACTS_VERBOSE("====================");
0203
0204       if (m_cfg.readCachedSurfaceInformation) {
0205         // add the surface information to the interaction this allows the
0206         // mapping to be speed up
0207         mInteraction.intersectionID = Acts::GeometryIdentifier((*m_sur_id)[is]);
0208         mInteraction.intersection =
0209             Acts::Vector3((*m_sur_x)[is], (*m_sur_y)[is], (*m_sur_z)[is]);
0210         mInteraction.pathCorrection = (*m_sur_pathCorrection)[is];
0211       } else {
0212         mInteraction.intersectionID = Acts::GeometryIdentifier();
0213         mInteraction.intersection = Acts::Vector3(0, 0, 0);
0214       }
0215       rmTrack.second.materialInteractions.push_back(std::move(mInteraction));
0216     }
0217     mtrackCollection[ib] = (std::move(rmTrack));
0218   }
0219
0220   // Write to the collection to the EventStore
0221   m_outputMaterialTracks(context, std::move(mtrackCollection));
0222
0223   // Return success flag
0224   return ProcessCode::SUCCESS;
0225 }
0226
0227 }  // namespace ActsExamples