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);
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     // necessary to guarantee that m_inputChain->GetV1() is valid for the
0103     // entire range
0104     m_inputChain->SetEstimate(m_events + 1);
0105
0106     m_entryNumbers.resize(m_events);
0107     m_inputChain->Draw("event_nr", "", "goff");
0108     RootUtility::stableSort(m_inputChain->GetEntries(), m_inputChain->GetV1(),
0109                             m_entryNumbers.data(), false);
0110   }
0111 }
0112
0113 std::pair<std::size_t, std::size_t> RootTrackSummaryReader::availableEvents()
0114     const {
0115   return {0u, m_events};
0116 }
0117
0118 RootTrackSummaryReader::~RootTrackSummaryReader() {
0119   delete m_multiTrajNr;
0120   delete m_subTrajNr;
0121   delete m_nStates;
0122   delete m_nMeasurements;
0123   delete m_nOutliers;
0124   delete m_nHoles;
0125   delete m_chi2Sum;
0126   delete m_NDF;
0127   delete m_measurementChi2;
0128   delete m_outlierChi2;
0129   delete m_measurementVolume;
0130   delete m_measurementLayer;
0131   delete m_outlierVolume;
0132   delete m_outlierLayer;
0133   delete m_majorityParticleId;
0134   delete m_nMajorityHits;
0135   delete m_t_charge;
0136   delete m_t_time;
0137   delete m_t_vx;
0138   delete m_t_vy;
0139   delete m_t_vz;
0140   delete m_t_px;
0141   delete m_t_py;
0142   delete m_t_pz;
0143   delete m_t_theta;
0144   delete m_t_phi;
0145   delete m_t_pT;
0146   delete m_t_eta;
0147   delete m_hasFittedParams;
0148   delete m_eLOC0_fit;
0149   delete m_eLOC1_fit;
0150   delete m_ePHI_fit;
0151   delete m_eTHETA_fit;
0152   delete m_eQOP_fit;
0153   delete m_eT_fit;
0154   delete m_err_eLOC0_fit;
0155   delete m_err_eLOC1_fit;
0156   delete m_err_ePHI_fit;
0157   delete m_err_eTHETA_fit;
0158   delete m_err_eQOP_fit;
0159   delete m_err_eT_fit;
0160 }
0161
0162 ProcessCode RootTrackSummaryReader::read(const AlgorithmContext& context) {
0163   ACTS_DEBUG("Trying to read recorded tracks.");
0164
0165   // read in the fitted track parameters and particles
0166   if (m_inputChain != nullptr && context.eventNumber < m_events) {
0167     // lock the mutex
0168     std::lock_guard<std::mutex> lock(m_read_mutex);
0169     // now read
0170
0171     std::shared_ptr<Acts::PerigeeSurface> perigeeSurface =
0172         Acts::Surface::makeShared<Acts::PerigeeSurface>(
0173             Acts::Vector3(0., 0., 0.));
0174
0175     // The collection to be written
0176     TrackParametersContainer trackParameterCollection;
0177     SimParticleContainer truthParticleCollection;
0178
0179     // Read the correct entry
0180     auto entry = m_entryNumbers.at(context.eventNumber);
0181     m_inputChain->GetEntry(entry);
0182     ACTS_INFO("Reading event: " << context.eventNumber
0183                                 << " stored as entry: " << entry);
0184
0185     unsigned int nTracks = m_eLOC0_fit->size();
0186     for (unsigned int i = 0; i < nTracks; i++) {
0187       Acts::BoundVector paramVec;
0188       paramVec << (*m_eLOC0_fit)[i], (*m_eLOC1_fit)[i], (*m_ePHI_fit)[i],
0189           (*m_eTHETA_fit)[i], (*m_eQOP_fit)[i], (*m_eT_fit)[i];
0190
0191       // Resolutions
0192       double resD0 = (*m_err_eLOC0_fit)[i];
0193       double resZ0 = (*m_err_eLOC1_fit)[i];
0194       double resPh = (*m_err_ePHI_fit)[i];
0195       double resTh = (*m_err_eTHETA_fit)[i];
0196       double resQp = (*m_err_eQOP_fit)[i];
0197       double resT = (*m_err_eT_fit)[i];
0198
0199       // Fill vector of track objects with simple covariance matrix
0200       Acts::BoundSquareMatrix covMat;
0201
0202       covMat << resD0 * resD0, 0., 0., 0., 0., 0., 0., resZ0 * resZ0, 0., 0.,
0203           0., 0., 0., 0., resPh * resPh, 0., 0., 0., 0., 0., 0., resTh * resTh,
0204           0., 0., 0., 0., 0., 0., resQp * resQp, 0., 0., 0., 0., 0., 0.,
0205           resT * resT;
0206
0207       // TODO we do not have a hypothesis at hand here. defaulting to pion
0208       trackParameterCollection.push_back(Acts::BoundTrackParameters(
0209           perigeeSurface, paramVec, std::move(covMat),
0210           Acts::ParticleHypothesis::pion()));
0211     }
0212
0213     unsigned int nTruthParticles = m_t_vx->size();
0214     for (unsigned int i = 0; i < nTruthParticles; i++) {
0215       SimParticleState truthParticle;
0216
0217       truthParticle.setPosition4((*m_t_vx)[i], (*m_t_vy)[i], (*m_t_vz)[i],
0218                                  (*m_t_time)[i]);
0219       truthParticle.setDirection((*m_t_px)[i], (*m_t_py)[i], (*m_t_pz)[i]);
0220       truthParticle.setParticleId((*m_majorityParticleId)[i]);
0221
0222       truthParticleCollection.insert(truthParticleCollection.end(),
0223                                      SimParticle(truthParticle, truthParticle));
0224     }
0225     // Write the collections to the EventStore
0226     m_outputTrackParameters(context, std::move(trackParameterCollection));
0227     m_outputParticles(context, std::move(truthParticleCollection));
0228   } else {
0229     ACTS_WARNING("Could not read in event.");
0230   }
0231   // Return success flag
0232   return ProcessCode::SUCCESS;
0233 }
0234
0235 }  // namespace ActsExamples