Root/src/RootMaterialTrackReader.cpp

0001 // This file is part of the Acts project.
0002 //
0003 // Copyright (C) 2017-2018 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 http://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 = new 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   std::cout << "The full chain has " << nentries << " entries for " << m_events
0087             << " events this corresponds to a batch size of: " << m_batchSize
0088             << std::endl;
0089
0090   // Sort the entry numbers of the events
0091   {
0092     m_entryNumbers.resize(nentries);
0093     m_inputChain->Draw("event_id", "", "goff");
0094     RootUtility::stableSort(m_inputChain->GetEntries(), m_inputChain->GetV1(),
0095                             m_entryNumbers.data(), false);
0096   }
0097
0098   m_outputMaterialTracks.initialize(m_cfg.outputMaterialTracks);
0099 }
0100
0101 RootMaterialTrackReader::~RootMaterialTrackReader() {
0102   delete m_inputChain;
0103
0104   delete m_step_x;
0105   delete m_step_y;
0106   delete m_step_z;
0107   delete m_step_dx;
0108   delete m_step_dy;
0109   delete m_step_dz;
0110   delete m_step_length;
0111   delete m_step_X0;
0112   delete m_step_L0;
0113   delete m_step_A;
0114   delete m_step_Z;
0115   delete m_step_rho;
0116
0117   delete m_sur_id;
0118   delete m_sur_x;
0119   delete m_sur_y;
0120   delete m_sur_z;
0121   delete m_sur_pathCorrection;
0122 }
0123
0124 std::string RootMaterialTrackReader::name() const {
0125   return "RootMaterialTrackReader";
0126 }
0127
0128 std::pair<std::size_t, std::size_t> RootMaterialTrackReader::availableEvents()
0129     const {
0130   return {0u, m_events};
0131 }
0132
0133 ProcessCode RootMaterialTrackReader::read(const AlgorithmContext& context) {
0134   ACTS_DEBUG("Trying to read recorded material from tracks.");
0135
0136   if (m_inputChain == nullptr || context.eventNumber >= m_events) {
0137     return ProcessCode::SUCCESS;
0138   }
0139
0140   // lock the mutex
0141   std::lock_guard<std::mutex> lock(m_read_mutex);
0142   // now read
0143
0144   // The collection to be written
0145   std::unordered_map<std::size_t, Acts::RecordedMaterialTrack> mtrackCollection;
0146
0147   // Loop over the entries for this event
0148   for (std::size_t ib = 0; ib < m_batchSize; ++ib) {
0149     // Read the correct entry: startEntry + ib
0150     auto entry = m_batchSize * context.eventNumber + ib;
0151     entry = m_entryNumbers.at(entry);
0152     ACTS_VERBOSE("Reading event: " << context.eventNumber
0153                                    << " with stored entry: " << entry);
0154     m_inputChain->GetEntry(entry);
0155
0156     Acts::RecordedMaterialTrack rmTrack;
0157     // Fill the position and momentum
0158     rmTrack.first.first = Acts::Vector3(m_v_x, m_v_y, m_v_z);
0159     rmTrack.first.second = Acts::Vector3(m_v_px, m_v_py, m_v_pz);
0160
0161     ACTS_VERBOSE("Track vertex:  " << rmTrack.first.first);
0162     ACTS_VERBOSE("Track momentum:" << rmTrack.first.second);
0163
0164     // Fill the individual steps
0165     std::size_t msteps = m_step_length->size();
0166     ACTS_VERBOSE("Reading " << msteps << " material steps.");
0167     rmTrack.second.materialInteractions.reserve(msteps);
0168     rmTrack.second.materialInX0 = 0.;
0169     rmTrack.second.materialInL0 = 0.;
0170
0171     for (std::size_t is = 0; is < msteps; ++is) {
0172       ACTS_VERBOSE("====================");
0173       ACTS_VERBOSE("[" << is + 1 << "/" << msteps << "] STEP INFORMATION: ");
0174
0175       double s = (*m_step_length)[is];
0176       if (s == 0) {
0177         ACTS_VERBOSE("invalid step length... skipping!");
0178         continue;
0179       }
0180
0181       double mX0 = (*m_step_X0)[is];
0182       double mL0 = (*m_step_L0)[is];
0183
0184       rmTrack.second.materialInX0 += s / mX0;
0185       rmTrack.second.materialInL0 += s / mL0;
0186       /// Fill the position & the material
0187       Acts::MaterialInteraction mInteraction;
0188       mInteraction.position =
0189           Acts::Vector3((*m_step_x)[is], (*m_step_y)[is], (*m_step_z)[is]);
0190       ACTS_VERBOSE("POSITION : " << (*m_step_x)[is] << ", " << (*m_step_y)[is]
0191                                  << ", " << (*m_step_z)[is]);
0192       mInteraction.direction =
0193           Acts::Vector3((*m_step_dx)[is], (*m_step_dy)[is], (*m_step_dz)[is]);
0194       ACTS_VERBOSE("DIRECTION: " << (*m_step_dx)[is] << ", " << (*m_step_dy)[is]
0195                                  << ", " << (*m_step_dz)[is]);
0196       mInteraction.materialSlab = Acts::MaterialSlab(
0197           Acts::Material::fromMassDensity(mX0, mL0, (*m_step_A)[is],
0198                                           (*m_step_Z)[is], (*m_step_rho)[is]),
0199           s);
0200       ACTS_VERBOSE("MATERIAL: " << mX0 << ", " << mL0 << ", " << (*m_step_A)[is]
0201                                 << ", " << (*m_step_Z)[is] << ", "
0202                                 << (*m_step_rho)[is]);
0203       ACTS_VERBOSE("====================");
0204
0205       if (m_cfg.readCachedSurfaceInformation) {
0206         // add the surface information to the interaction this allows the
0207         // mapping to be speed up
0208         mInteraction.intersectionID = Acts::GeometryIdentifier((*m_sur_id)[is]);
0209         mInteraction.intersection =
0210             Acts::Vector3((*m_sur_x)[is], (*m_sur_y)[is], (*m_sur_z)[is]);
0211         mInteraction.pathCorrection = (*m_sur_pathCorrection)[is];
0212       } else {
0213         mInteraction.intersectionID = Acts::GeometryIdentifier();
0214         mInteraction.intersection = Acts::Vector3(0, 0, 0);
0215       }
0216       rmTrack.second.materialInteractions.push_back(std::move(mInteraction));
0217     }
0218     mtrackCollection[ib] = (std::move(rmTrack));
0219   }
0220
0221   // Write to the collection to the EventStore
0222   m_outputMaterialTracks(context, std::move(mtrackCollection));
0223
0224   // Return success flag
0225   return ProcessCode::SUCCESS;
0226 }
0227
0228 }  // namespace ActsExamples