Root/src/RootTrackSummaryReader.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/RootTrackSummaryReader.hpp"
0010
0011 #include "Acts/Definitions/Algebra.hpp"
0012 #include "Acts/Definitions/TrackParametrization.hpp"
0013 #include "Acts/EventData/TrackParameters.hpp"
0014 #include "Acts/Surfaces/PerigeeSurface.hpp"
0015 #include "Acts/Surfaces/Surface.hpp"
0016 #include "Acts/Utilities/Logger.hpp"
0017 #include "ActsExamples/EventData/SimParticle.hpp"
0018 #include "ActsExamples/Framework/AlgorithmContext.hpp"
0019 #include "ActsExamples/Io/Root/RootUtility.hpp"
0020 #include "ActsFatras/EventData/Particle.hpp"
0021
0022 #include <iostream>
0023 #include <stdexcept>
0024
0025 #include <TChain.h>
0026
0027 namespace ActsExamples {
0028
0029 RootTrackSummaryReader::RootTrackSummaryReader(
0030     const RootTrackSummaryReader::Config& config, Acts::Logging::Level level)
0031     : IReader(),
0032       m_logger{Acts::getDefaultLogger(name(), level)},
0033       m_cfg(config) {
0034   m_inputChain = std::make_unique<TChain>(m_cfg.treeName.c_str());
0035
0036   if (m_cfg.filePath.empty()) {
0037     throw std::invalid_argument("Missing input filename");
0038   }
0039
0040   m_outputTrackParameters.initialize(m_cfg.outputTracks);
0041   m_outputParticles.initialize(m_cfg.outputParticles);
0042
0043   // Set the branches
0044   m_inputChain->SetBranchAddress("event_nr", &m_eventNr);
0045   m_inputChain->SetBranchAddress("multiTraj_nr", &m_multiTrajNr.get());
0046   m_inputChain->SetBranchAddress("subTraj_nr", &m_subTrajNr.get());
0047
0048   // These info is not really stored in the event store, but still read in
0049   m_inputChain->SetBranchAddress("nStates", &m_nStates.get());
0050   m_inputChain->SetBranchAddress("nMeasurements", &m_nMeasurements.get());
0051   m_inputChain->SetBranchAddress("nOutliers", &m_nOutliers.get());
0052   m_inputChain->SetBranchAddress("nHoles", &m_nHoles.get());
0053   m_inputChain->SetBranchAddress("chi2Sum", &m_chi2Sum.get());
0054   m_inputChain->SetBranchAddress("NDF", &m_NDF.get());
0055   m_inputChain->SetBranchAddress("measurementChi2", &m_measurementChi2.get());
0056   m_inputChain->SetBranchAddress("outlierChi2", &m_outlierChi2.get());
0057   m_inputChain->SetBranchAddress("measurementVolume",
0058                                  &m_measurementVolume.get());
0059   m_inputChain->SetBranchAddress("measurementLayer", &m_measurementLayer.get());
0060   m_inputChain->SetBranchAddress("outlierVolume", &m_outlierVolume.get());
0061   m_inputChain->SetBranchAddress("outlierLayer", &m_outlierLayer.get());
0062
0063   if (m_inputChain->GetBranch("majorityParticleId") != nullptr) {
0064     m_hasCombinedMajorityParticleId = true;
0065     m_majorityParticleId.allocate();
0066     m_inputChain->SetBranchAddress("majorityParticleId",
0067                                    &m_majorityParticleId.get());
0068   } else {
0069     m_hasCombinedMajorityParticleId = false;
0070     m_majorityParticleVertexPrimary.allocate();
0071     m_majorityParticleVertexSecondary.allocate();
0072     m_majorityParticleParticle.allocate();
0073     m_majorityParticleGeneration.allocate();
0074     m_majorityParticleSubParticle.allocate();
0075     m_inputChain->SetBranchAddress("majorityParticleId_vertex_primary",
0076                                    &m_majorityParticleVertexPrimary.get());
0077     m_inputChain->SetBranchAddress("majorityParticleId_vertex_secondary",
0078                                    &m_majorityParticleVertexSecondary.get());
0079     m_inputChain->SetBranchAddress("majorityParticleId_particle",
0080                                    &m_majorityParticleParticle.get());
0081     m_inputChain->SetBranchAddress("majorityParticleId_generation",
0082                                    &m_majorityParticleGeneration.get());
0083     m_inputChain->SetBranchAddress("majorityParticleId_sub_particle",
0084                                    &m_majorityParticleSubParticle.get());
0085   }
0086   m_inputChain->SetBranchAddress("nMajorityHits", &m_nMajorityHits.get());
0087   m_inputChain->SetBranchAddress("t_charge", &m_t_charge.get());
0088   m_inputChain->SetBranchAddress("t_time", &m_t_time.get());
0089   m_inputChain->SetBranchAddress("t_vx", &m_t_vx.get());
0090   m_inputChain->SetBranchAddress("t_vy", &m_t_vy.get());
0091   m_inputChain->SetBranchAddress("t_vz", &m_t_vz.get());
0092   m_inputChain->SetBranchAddress("t_px", &m_t_px.get());
0093   m_inputChain->SetBranchAddress("t_py", &m_t_py.get());
0094   m_inputChain->SetBranchAddress("t_pz", &m_t_pz.get());
0095   m_inputChain->SetBranchAddress("t_theta", &m_t_theta.get());
0096   m_inputChain->SetBranchAddress("t_phi", &m_t_phi.get());
0097   m_inputChain->SetBranchAddress("t_eta", &m_t_eta.get());
0098   m_inputChain->SetBranchAddress("t_pT", &m_t_pT.get());
0099
0100   m_inputChain->SetBranchAddress("hasFittedParams", &m_hasFittedParams.get());
0101   m_inputChain->SetBranchAddress("eLOC0_fit", &m_eLOC0_fit.get());
0102   m_inputChain->SetBranchAddress("eLOC1_fit", &m_eLOC1_fit.get());
0103   m_inputChain->SetBranchAddress("ePHI_fit", &m_ePHI_fit.get());
0104   m_inputChain->SetBranchAddress("eTHETA_fit", &m_eTHETA_fit.get());
0105   m_inputChain->SetBranchAddress("eQOP_fit", &m_eQOP_fit.get());
0106   m_inputChain->SetBranchAddress("eT_fit", &m_eT_fit.get());
0107   m_inputChain->SetBranchAddress("err_eLOC0_fit", &m_err_eLOC0_fit.get());
0108   m_inputChain->SetBranchAddress("err_eLOC1_fit", &m_err_eLOC1_fit.get());
0109   m_inputChain->SetBranchAddress("err_ePHI_fit", &m_err_ePHI_fit.get());
0110   m_inputChain->SetBranchAddress("err_eTHETA_fit", &m_err_eTHETA_fit.get());
0111   m_inputChain->SetBranchAddress("err_eQOP_fit", &m_err_eQOP_fit.get());
0112   m_inputChain->SetBranchAddress("err_eT_fit", &m_err_eT_fit.get());
0113
0114   auto path = m_cfg.filePath;
0115
0116   // add file to the input chain
0117   m_inputChain->Add(path.c_str());
0118   ACTS_DEBUG("Adding File " << path << " to tree '" << m_cfg.treeName << "'.");
0119
0120   m_events = m_inputChain->GetEntries();
0121   ACTS_DEBUG("The full chain has " << m_events << " entries.");
0122
0123   // Sort the entry numbers of the events
0124   {
0125     // necessary to guarantee that m_inputChain->GetV1() is valid for the
0126     // entire range
0127     m_inputChain->SetEstimate(m_events + 1);
0128
0129     m_entryNumbers.resize(m_events);
0130     m_inputChain->Draw("event_nr", "", "goff");
0131     RootUtility::stableSort(m_inputChain->GetEntries(), m_inputChain->GetV1(),
0132                             m_entryNumbers.data(), false);
0133   }
0134 }
0135
0136 std::pair<std::size_t, std::size_t> RootTrackSummaryReader::availableEvents()
0137     const {
0138   return {0u, m_events};
0139 }
0140
0141 ProcessCode RootTrackSummaryReader::read(const AlgorithmContext& context) {
0142   ACTS_DEBUG("Trying to read recorded tracks.");
0143
0144   // read in the fitted track parameters and particles
0145   if (m_inputChain != nullptr && context.eventNumber < m_events) {
0146     // lock the mutex
0147     std::lock_guard<std::mutex> lock(m_read_mutex);
0148     // now read
0149
0150     std::shared_ptr<Acts::PerigeeSurface> perigeeSurface =
0151         Acts::Surface::makeShared<Acts::PerigeeSurface>(
0152             Acts::Vector3(0., 0., 0.));
0153
0154     // The collection to be written
0155     TrackParametersContainer trackParameterCollection;
0156     SimParticleContainer truthParticleCollection;
0157
0158     // Read the correct entry
0159     auto entry = m_entryNumbers.at(context.eventNumber);
0160     m_inputChain->GetEntry(entry);
0161     ACTS_INFO("Reading event: " << context.eventNumber
0162                                 << " stored as entry: " << entry);
0163
0164     unsigned int nTracks = m_eLOC0_fit->size();
0165     for (unsigned int i = 0; i < nTracks; i++) {
0166       Acts::BoundVector paramVec;
0167       paramVec << (*m_eLOC0_fit)[i], (*m_eLOC1_fit)[i], (*m_ePHI_fit)[i],
0168           (*m_eTHETA_fit)[i], (*m_eQOP_fit)[i], (*m_eT_fit)[i];
0169
0170       // Resolutions
0171       double resD0 = (*m_err_eLOC0_fit)[i];
0172       double resZ0 = (*m_err_eLOC1_fit)[i];
0173       double resPh = (*m_err_ePHI_fit)[i];
0174       double resTh = (*m_err_eTHETA_fit)[i];
0175       double resQp = (*m_err_eQOP_fit)[i];
0176       double resT = (*m_err_eT_fit)[i];
0177
0178       // Fill vector of track objects with simple covariance matrix
0179       Acts::BoundSquareMatrix covMat;
0180
0181       covMat << resD0 * resD0, 0., 0., 0., 0., 0., 0., resZ0 * resZ0, 0., 0.,
0182           0., 0., 0., 0., resPh * resPh, 0., 0., 0., 0., 0., 0., resTh * resTh,
0183           0., 0., 0., 0., 0., 0., resQp * resQp, 0., 0., 0., 0., 0., 0.,
0184           resT * resT;
0185
0186       // TODO we do not have a hypothesis at hand here. defaulting to pion
0187       trackParameterCollection.push_back(Acts::BoundTrackParameters(
0188           perigeeSurface, paramVec, std::move(covMat),
0189           Acts::ParticleHypothesis::pion()));
0190     }
0191
0192     unsigned int nTruthParticles = m_t_vx->size();
0193     for (unsigned int i = 0; i < nTruthParticles; i++) {
0194       SimParticleState truthParticle;
0195
0196       truthParticle.setPosition4((*m_t_vx)[i], (*m_t_vy)[i], (*m_t_vz)[i],
0197                                  (*m_t_time)[i]);
0198       truthParticle.setDirection((*m_t_px)[i], (*m_t_py)[i], (*m_t_pz)[i]);
0199       auto makeBarcode = [](std::uint32_t vp, std::uint32_t vs,
0200                             std::uint32_t particle, std::uint32_t generation,
0201                             std::uint32_t subParticle) {
0202         SimBarcode barcode = SimBarcode::Invalid();
0203         return barcode
0204             .withVertexPrimary(static_cast<SimBarcode::PrimaryVertexId>(vp))
0205             .withVertexSecondary(static_cast<SimBarcode::SecondaryVertexId>(vs))
0206             .withParticle(static_cast<SimBarcode::ParticleId>(particle))
0207             .withGeneration(static_cast<SimBarcode::GenerationId>(generation))
0208             .withSubParticle(
0209                 static_cast<SimBarcode::SubParticleId>(subParticle));
0210       };
0211
0212       SimBarcode barcode = SimBarcode::Invalid();
0213       if (m_hasCombinedMajorityParticleId && m_majorityParticleId.hasValue()) {
0214         const auto& components = (*m_majorityParticleId)[i];
0215         auto comp = [&](std::size_t idx) -> std::uint32_t {
0216           return (components.size() > idx) ? components[idx] : 0u;
0217         };
0218         barcode = makeBarcode(comp(0), comp(1), comp(2), comp(3), comp(4));
0219       } else {
0220         auto safeAt = [](const auto& branch, std::size_t idx) -> std::uint32_t {
0221           const auto* vec = branch.get();
0222           return (vec != nullptr && vec->size() > idx) ? vec->at(idx) : 0u;
0223         };
0224         barcode = makeBarcode(safeAt(m_majorityParticleVertexPrimary, i),
0225                               safeAt(m_majorityParticleVertexSecondary, i),
0226                               safeAt(m_majorityParticleParticle, i),
0227                               safeAt(m_majorityParticleGeneration, i),
0228                               safeAt(m_majorityParticleSubParticle, i));
0229       }
0230       truthParticle.setParticleId(barcode);
0231
0232       truthParticleCollection.insert(truthParticleCollection.end(),
0233                                      SimParticle(truthParticle, truthParticle));
0234     }
0235     // Write the collections to the EventStore
0236     m_outputTrackParameters(context, std::move(trackParameterCollection));
0237     m_outputParticles(context, std::move(truthParticleCollection));
0238   } else {
0239     ACTS_WARNING("Could not read in event.");
0240   }
0241   // Return success flag
0242   return ProcessCode::SUCCESS;
0243 }
0244
0245 }  // namespace ActsExamples