Root/src/RootTrackSummaryReader.cpp

0001 // This file is part of the Acts project.
0002 //
0003 // Copyright (C) 2021 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/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 = new 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);
0046   m_inputChain->SetBranchAddress("subTraj_nr", &m_subTrajNr);
0047
0048   // These info is not really stored in the event store, but still read in
0049   m_inputChain->SetBranchAddress("nStates", &m_nStates);
0050   m_inputChain->SetBranchAddress("nMeasurements", &m_nMeasurements);
0051   m_inputChain->SetBranchAddress("nOutliers", &m_nOutliers);
0052   m_inputChain->SetBranchAddress("nHoles", &m_nHoles);
0053   m_inputChain->SetBranchAddress("chi2Sum", &m_chi2Sum);
0054   m_inputChain->SetBranchAddress("NDF", &m_NDF);
0055   m_inputChain->SetBranchAddress("measurementChi2", &m_measurementChi2);
0056   m_inputChain->SetBranchAddress("outlierChi2", &m_outlierChi2);
0057   m_inputChain->SetBranchAddress("measurementVolume", &m_measurementVolume);
0058   m_inputChain->SetBranchAddress("measurementLayer", &m_measurementLayer);
0059   m_inputChain->SetBranchAddress("outlierVolume", &m_outlierVolume);
0060   m_inputChain->SetBranchAddress("outlierLayer", &m_outlierLayer);
0061
0062   m_inputChain->SetBranchAddress("majorityParticleId", &m_majorityParticleId);
0063   m_inputChain->SetBranchAddress("nMajorityHits", &m_nMajorityHits);
0064   m_inputChain->SetBranchAddress("t_charge", &m_t_charge);
0065   m_inputChain->SetBranchAddress("t_time", &m_t_time);
0066   m_inputChain->SetBranchAddress("t_vx", &m_t_vx);
0067   m_inputChain->SetBranchAddress("t_vy", &m_t_vy);
0068   m_inputChain->SetBranchAddress("t_vz", &m_t_vz);
0069   m_inputChain->SetBranchAddress("t_px", &m_t_px);
0070   m_inputChain->SetBranchAddress("t_py", &m_t_py);
0071   m_inputChain->SetBranchAddress("t_pz", &m_t_pz);
0072   m_inputChain->SetBranchAddress("t_theta", &m_t_theta);
0073   m_inputChain->SetBranchAddress("t_phi", &m_t_phi);
0074   m_inputChain->SetBranchAddress("t_eta", &m_t_eta);
0075   m_inputChain->SetBranchAddress("t_pT", &m_t_pT);
0076
0077   m_inputChain->SetBranchAddress("hasFittedParams", &m_hasFittedParams);
0078   m_inputChain->SetBranchAddress("eLOC0_fit", &m_eLOC0_fit);
0079   m_inputChain->SetBranchAddress("eLOC1_fit", &m_eLOC1_fit);
0080   m_inputChain->SetBranchAddress("ePHI_fit", &m_ePHI_fit);
0081   m_inputChain->SetBranchAddress("eTHETA_fit", &m_eTHETA_fit);
0082   m_inputChain->SetBranchAddress("eQOP_fit", &m_eQOP_fit);
0083   m_inputChain->SetBranchAddress("eT_fit", &m_eT_fit);
0084   m_inputChain->SetBranchAddress("err_eLOC0_fit", &m_err_eLOC0_fit);
0085   m_inputChain->SetBranchAddress("err_eLOC1_fit", &m_err_eLOC1_fit);
0086   m_inputChain->SetBranchAddress("err_ePHI_fit", &m_err_ePHI_fit);
0087   m_inputChain->SetBranchAddress("err_eTHETA_fit", &m_err_eTHETA_fit);
0088   m_inputChain->SetBranchAddress("err_eQOP_fit", &m_err_eQOP_fit);
0089   m_inputChain->SetBranchAddress("err_eT_fit", &m_err_eT_fit);
0090
0091   auto path = m_cfg.filePath;
0092
0093   // add file to the input chain
0094   m_inputChain->Add(path.c_str());
0095   ACTS_DEBUG("Adding File " << path << " to tree '" << m_cfg.treeName << "'.");
0096
0097   m_events = m_inputChain->GetEntries();
0098   ACTS_DEBUG("The full chain has " << m_events << " entries.");
0099
0100   // Sort the entry numbers of the events
0101   {
0102     m_entryNumbers.resize(m_events);
0103     m_inputChain->Draw("event_nr", "", "goff");
0104     RootUtility::stableSort(m_inputChain->GetEntries(), m_inputChain->GetV1(),
0105                             m_entryNumbers.data(), false);
0106   }
0107 }
0108
0109 std::pair<std::size_t, std::size_t> RootTrackSummaryReader::availableEvents()
0110     const {
0111   return {0u, m_events};
0112 }
0113
0114 RootTrackSummaryReader::~RootTrackSummaryReader() {
0115   delete m_multiTrajNr;
0116   delete m_subTrajNr;
0117   delete m_nStates;
0118   delete m_nMeasurements;
0119   delete m_nOutliers;
0120   delete m_nHoles;
0121   delete m_chi2Sum;
0122   delete m_NDF;
0123   delete m_measurementChi2;
0124   delete m_outlierChi2;
0125   delete m_measurementVolume;
0126   delete m_measurementLayer;
0127   delete m_outlierVolume;
0128   delete m_outlierLayer;
0129   delete m_majorityParticleId;
0130   delete m_nMajorityHits;
0131   delete m_t_charge;
0132   delete m_t_time;
0133   delete m_t_vx;
0134   delete m_t_vy;
0135   delete m_t_vz;
0136   delete m_t_px;
0137   delete m_t_py;
0138   delete m_t_pz;
0139   delete m_t_theta;
0140   delete m_t_phi;
0141   delete m_t_pT;
0142   delete m_t_eta;
0143   delete m_hasFittedParams;
0144   delete m_eLOC0_fit;
0145   delete m_eLOC1_fit;
0146   delete m_ePHI_fit;
0147   delete m_eTHETA_fit;
0148   delete m_eQOP_fit;
0149   delete m_eT_fit;
0150   delete m_err_eLOC0_fit;
0151   delete m_err_eLOC1_fit;
0152   delete m_err_ePHI_fit;
0153   delete m_err_eTHETA_fit;
0154   delete m_err_eQOP_fit;
0155   delete m_err_eT_fit;
0156 }
0157
0158 ProcessCode RootTrackSummaryReader::read(const AlgorithmContext& context) {
0159   ACTS_DEBUG("Trying to read recorded tracks.");
0160
0161   // read in the fitted track parameters and particles
0162   if (m_inputChain != nullptr && context.eventNumber < m_events) {
0163     // lock the mutex
0164     std::lock_guard<std::mutex> lock(m_read_mutex);
0165     // now read
0166
0167     std::shared_ptr<Acts::PerigeeSurface> perigeeSurface =
0168         Acts::Surface::makeShared<Acts::PerigeeSurface>(
0169             Acts::Vector3(0., 0., 0.));
0170
0171     // The collection to be written
0172     TrackParametersContainer trackParameterCollection;
0173     SimParticleContainer truthParticleCollection;
0174
0175     // Read the correct entry
0176     auto entry = m_entryNumbers.at(context.eventNumber);
0177     m_inputChain->GetEntry(entry);
0178     ACTS_INFO("Reading event: " << context.eventNumber
0179                                 << " stored as entry: " << entry);
0180
0181     unsigned int nTracks = m_eLOC0_fit->size();
0182     for (unsigned int i = 0; i < nTracks; i++) {
0183       Acts::BoundVector paramVec;
0184       paramVec << (*m_eLOC0_fit)[i], (*m_eLOC1_fit)[i], (*m_ePHI_fit)[i],
0185           (*m_eTHETA_fit)[i], (*m_eQOP_fit)[i], (*m_eT_fit)[i];
0186
0187       // Resolutions
0188       double resD0 = (*m_err_eLOC0_fit)[i];
0189       double resZ0 = (*m_err_eLOC1_fit)[i];
0190       double resPh = (*m_err_ePHI_fit)[i];
0191       double resTh = (*m_err_eTHETA_fit)[i];
0192       double resQp = (*m_err_eQOP_fit)[i];
0193       double resT = (*m_err_eT_fit)[i];
0194
0195       // Fill vector of track objects with simple covariance matrix
0196       Acts::BoundSquareMatrix covMat;
0197
0198       covMat << resD0 * resD0, 0., 0., 0., 0., 0., 0., resZ0 * resZ0, 0., 0.,
0199           0., 0., 0., 0., resPh * resPh, 0., 0., 0., 0., 0., 0., resTh * resTh,
0200           0., 0., 0., 0., 0., 0., resQp * resQp, 0., 0., 0., 0., 0., 0.,
0201           resT * resT;
0202
0203       // TODO we do not have a hypothesis at hand here. defaulting to pion
0204       trackParameterCollection.push_back(Acts::BoundTrackParameters(
0205           perigeeSurface, paramVec, std::move(covMat),
0206           Acts::ParticleHypothesis::pion()));
0207     }
0208
0209     unsigned int nTruthParticles = m_t_vx->size();
0210     for (unsigned int i = 0; i < nTruthParticles; i++) {
0211       ActsFatras::Particle truthParticle;
0212
0213       truthParticle.setPosition4((*m_t_vx)[i], (*m_t_vy)[i], (*m_t_vz)[i],
0214                                  (*m_t_time)[i]);
0215       truthParticle.setDirection((*m_t_px)[i], (*m_t_py)[i], (*m_t_pz)[i]);
0216       truthParticle.setParticleId((*m_majorityParticleId)[i]);
0217
0218       truthParticleCollection.insert(truthParticleCollection.end(),
0219                                      truthParticle);
0220     }
0221     // Write the collections to the EventStore
0222     m_outputTrackParameters(context, std::move(trackParameterCollection));
0223     m_outputParticles(context, std::move(truthParticleCollection));
0224   } else {
0225     ACTS_WARNING("Could not read in event.");
0226   }
0227   // Return success flag
0228   return ProcessCode::SUCCESS;
0229 }
0230
0231 }  // namespace ActsExamples