Root/src/RootAthenaDumpReader.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/RootAthenaDumpReader.hpp"
0010
0011 #include "Acts/Definitions/Units.hpp"
0012 #include "Acts/EventData/SourceLink.hpp"
0013 #include "Acts/Geometry/GeometryIdentifier.hpp"
0014 #include "Acts/Geometry/TrackingGeometry.hpp"
0015 #include "Acts/Utilities/Zip.hpp"
0016 #include "ActsExamples/Digitization/MeasurementCreation.hpp"
0017 #include "ActsExamples/EventData/Cluster.hpp"
0018 #include "ActsExamples/EventData/IndexSourceLink.hpp"
0019 #include "ActsExamples/EventData/SimParticle.hpp"
0020 #include "ActsExamples/EventData/SpacePoint.hpp"
0021
0022 #include <algorithm>
0023 #include <numeric>
0024
0025 #include <TChain.h>
0026 #include <boost/container/static_vector.hpp>
0027
0028 using namespace Acts::UnitLiterals;
0029
0030 namespace ActsExamples {
0031
0032 namespace {
0033
0034 /// In cases when there is built up a particle collection in an iterative way it
0035 /// can be way faster to build up a vector and afterwards use a special
0036 /// constructor to speed up the set creation.
0037 inline auto particleVectorToSet(std::vector<SimParticle>& particles) {
0038   using namespace ActsExamples;
0039   auto cmp = [](const auto& a, const auto& b) {
0040     return a.particleId() == b.particleId();
0041   };
0042
0043   std::ranges::sort(particles, detail::CompareParticleId{});
0044   particles.erase(std::unique(particles.begin(), particles.end(), cmp),
0045                   particles.end());
0046
0047   return SimParticleContainer(boost::container::ordered_unique_range_t{},
0048                               particles.begin(), particles.end());
0049 }
0050
0051 }  // namespace
0052
0053 enum SpacePointType { ePixel = 1, eStrip = 2 };
0054
0055 RootAthenaDumpReader::RootAthenaDumpReader(
0056     const RootAthenaDumpReader::Config& config, Acts::Logging::Level level)
0057     : IReader(),
0058       m_cfg(config),
0059       m_logger(Acts::getDefaultLogger(name(), level)) {
0060   if (m_cfg.inputfiles.empty()) {
0061     throw std::invalid_argument("Empty input file list");
0062   }
0063   if (m_cfg.treename.empty()) {
0064     throw std::invalid_argument("Missing tree name");
0065   }
0066
0067   m_inputchain = std::make_shared<TChain>(m_cfg.treename.c_str());
0068
0069   m_outputPixelSpacePoints.initialize(m_cfg.outputPixelSpacePoints);
0070   m_outputStripSpacePoints.initialize(m_cfg.outputStripSpacePoints);
0071   m_outputSpacePoints.initialize(m_cfg.outputSpacePoints);
0072   if (!m_cfg.onlySpacePoints) {
0073     m_outputMeasurements.initialize(m_cfg.outputMeasurements);
0074     m_outputMeasurementSubset.initialize(m_cfg.outputMeasurementSubset);
0075     m_outputClusters.initialize(m_cfg.outputClusters);
0076     if (!m_cfg.noTruth) {
0077       m_outputParticles.initialize(m_cfg.outputParticles);
0078       m_outputMeasParticleMap.initialize(m_cfg.outputMeasurementParticlesMap);
0079       m_outputParticleMeasMap.initialize(m_cfg.outputParticleMeasurementsMap);
0080     }
0081   }
0082
0083   if (m_inputchain->GetBranch("SPtopStripDirection") == nullptr) {
0084     ACTS_INFO("Additional SP strip features not available");
0085     m_haveStripFeatures = false;
0086   }
0087
0088   if (m_inputchain->GetBranch("CLetas") == nullptr) {
0089     ACTS_INFO("Cluster cell data (CLetas/CLphis/CLtots) not available");
0090     m_haveCellData = false;
0091   }
0092
0093   // Set the branches
0094   m_inputchain->SetBranchAddress("run_number", &run_number);
0095   m_inputchain->SetBranchAddress("event_number", &event_number);
0096   m_inputchain->SetBranchAddress("nSE", &nSE);
0097   m_inputchain->SetBranchAddress("SEID", SEID);
0098
0099   // Cluster features
0100   m_inputchain->SetBranchAddress("nCL", &nCL);
0101   m_inputchain->SetBranchAddress("CLindex", CLindex);
0102   m_inputchain->SetBranchAddress("CLhardware", &CLhardware.get());
0103   m_inputchain->SetBranchAddress("CLx", CLx);
0104   m_inputchain->SetBranchAddress("CLy", CLy);
0105   m_inputchain->SetBranchAddress("CLz", CLz);
0106   m_inputchain->SetBranchAddress("CLbarrel_endcap", CLbarrel_endcap);
0107   m_inputchain->SetBranchAddress("CLlayer_disk", CLlayer_disk);
0108   m_inputchain->SetBranchAddress("CLeta_module", CLeta_module);
0109   m_inputchain->SetBranchAddress("CLphi_module", CLphi_module);
0110   m_inputchain->SetBranchAddress("CLside", CLside);
0111   m_inputchain->SetBranchAddress("CLmoduleID", CLmoduleID);
0112   if (m_haveCellData) {
0113     m_inputchain->SetBranchAddress("CLphis", &CLphis.get());
0114     m_inputchain->SetBranchAddress("CLetas", &CLetas.get());
0115     m_inputchain->SetBranchAddress("CLtots", &CLtots.get());
0116   }
0117   m_inputchain->SetBranchAddress("CLloc_direction1", CLloc_direction1);
0118   m_inputchain->SetBranchAddress("CLloc_direction2", CLloc_direction2);
0119   m_inputchain->SetBranchAddress("CLloc_direction3", CLloc_direction3);
0120   m_inputchain->SetBranchAddress("CLJan_loc_direction1", CLJan_loc_direction1);
0121   m_inputchain->SetBranchAddress("CLJan_loc_direction2", CLJan_loc_direction2);
0122   m_inputchain->SetBranchAddress("CLJan_loc_direction3", CLJan_loc_direction3);
0123   m_inputchain->SetBranchAddress("CLpixel_count", CLpixel_count);
0124   m_inputchain->SetBranchAddress("CLcharge_count", CLcharge_count);
0125   m_inputchain->SetBranchAddress("CLloc_eta", CLloc_eta);
0126   m_inputchain->SetBranchAddress("CLloc_phi", CLloc_phi);
0127   m_inputchain->SetBranchAddress("CLglob_eta", CLglob_eta);
0128   m_inputchain->SetBranchAddress("CLglob_phi", CLglob_phi);
0129   m_inputchain->SetBranchAddress("CLeta_angle", CLeta_angle);
0130   m_inputchain->SetBranchAddress("CLphi_angle", CLphi_angle);
0131   m_inputchain->SetBranchAddress("CLnorm_x", CLnorm_x);
0132   m_inputchain->SetBranchAddress("CLnorm_y", CLnorm_y);
0133   m_inputchain->SetBranchAddress("CLnorm_z", CLnorm_z);
0134   m_inputchain->SetBranchAddress("CLlocal_cov", &CLlocal_cov.get());
0135   if (!m_cfg.noTruth) {
0136     m_inputchain->SetBranchAddress("CLparticleLink_eventIndex",
0137                                    &CLparticleLink_eventIndex.get());
0138     m_inputchain->SetBranchAddress("CLparticleLink_barcode",
0139                                    &CLparticleLink_barcode.get());
0140     m_inputchain->SetBranchAddress("CLbarcodesLinked", &CLbarcodesLinked.get());
0141     m_inputchain->SetBranchAddress("CLparticle_charge",
0142                                    &CLparticle_charge.get());
0143   }
0144
0145   // Particle features
0146   if (!m_cfg.noTruth) {
0147     m_inputchain->SetBranchAddress("nPartEVT", &nPartEVT);
0148     m_inputchain->SetBranchAddress("Part_event_number", Part_event_number);
0149     m_inputchain->SetBranchAddress("Part_barcode", Part_barcode);
0150     m_inputchain->SetBranchAddress("Part_px", Part_px);
0151     m_inputchain->SetBranchAddress("Part_py", Part_py);
0152     m_inputchain->SetBranchAddress("Part_pz", Part_pz);
0153     m_inputchain->SetBranchAddress("Part_pt", Part_pt);
0154     m_inputchain->SetBranchAddress("Part_eta", Part_eta);
0155     m_inputchain->SetBranchAddress("Part_vx", Part_vx);
0156     m_inputchain->SetBranchAddress("Part_vy", Part_vy);
0157     m_inputchain->SetBranchAddress("Part_vz", Part_vz);
0158     m_inputchain->SetBranchAddress("Part_radius", Part_radius);
0159     m_inputchain->SetBranchAddress("Part_status", Part_status);
0160     m_inputchain->SetBranchAddress("Part_charge", Part_charge);
0161     m_inputchain->SetBranchAddress("Part_pdg_id", Part_pdg_id);
0162     m_inputchain->SetBranchAddress("Part_passed", Part_passed);
0163     m_inputchain->SetBranchAddress("Part_vProdNin", Part_vProdNin);
0164     m_inputchain->SetBranchAddress("Part_vProdNout", Part_vProdNout);
0165     m_inputchain->SetBranchAddress("Part_vProdStatus", Part_vProdStatus);
0166     m_inputchain->SetBranchAddress("Part_vProdBarcode", Part_vProdBarcode);
0167     m_inputchain->SetBranchAddress("Part_vParentID", &Part_vParentID.get());
0168     m_inputchain->SetBranchAddress("Part_vParentBarcode",
0169                                    &Part_vParentBarcode.get());
0170   }
0171
0172   // Space point features
0173   m_inputchain->SetBranchAddress("nSP", &nSP);
0174   m_inputchain->SetBranchAddress("SPindex", SPindex);
0175   m_inputchain->SetBranchAddress("SPx", SPx);
0176   m_inputchain->SetBranchAddress("SPy", SPy);
0177   m_inputchain->SetBranchAddress("SPz", SPz);
0178   m_inputchain->SetBranchAddress("SPCL1_index", SPCL1_index);
0179   m_inputchain->SetBranchAddress("SPCL2_index", SPCL2_index);
0180   m_inputchain->SetBranchAddress("SPisOverlap", SPisOverlap);
0181   if (m_haveStripFeatures) {
0182     m_inputchain->SetBranchAddress("SPradius", SPradius);
0183     m_inputchain->SetBranchAddress("SPcovr", SPcovr);
0184     m_inputchain->SetBranchAddress("SPcovz", SPcovz);
0185     m_inputchain->SetBranchAddress("SPhl_topstrip", SPhl_topstrip);
0186     m_inputchain->SetBranchAddress("SPhl_botstrip", SPhl_botstrip);
0187     m_inputchain->SetBranchAddress("SPtopStripDirection",
0188                                    &SPtopStripDirection.get());
0189     m_inputchain->SetBranchAddress("SPbottomStripDirection",
0190                                    &SPbottomStripDirection.get());
0191     m_inputchain->SetBranchAddress("SPstripCenterDistance",
0192                                    &SPstripCenterDistance.get());
0193     m_inputchain->SetBranchAddress("SPtopStripCenterPosition",
0194                                    &SPtopStripCenterPosition.get());
0195   }
0196
0197   // These quantities are not used currently and thus commented out
0198   // I would like to keep the code, since it is always a pain to write it
0199   /*
0200   m_inputchain->SetBranchAddress("nTRK", &nTRK);
0201   m_inputchain->SetBranchAddress("TRKindex", TRKindex);
0202   m_inputchain->SetBranchAddress("TRKtrack_fitter", TRKtrack_fitter);
0203   m_inputchain->SetBranchAddress("TRKparticle_hypothesis",
0204                                  TRKparticle_hypothesis);
0205   m_inputchain->SetBranchAddress("TRKproperties", &TRKproperties);
0206   m_inputchain->SetBranchAddress("TRKpattern", &TRKpattern);
0207   m_inputchain->SetBranchAddress("TRKndof", TRKndof);
0208   m_inputchain->SetBranchAddress("TRKmot", TRKmot);
0209   m_inputchain->SetBranchAddress("TRKoot", TRKoot);
0210   m_inputchain->SetBranchAddress("TRKchiSq", TRKchiSq);
0211   m_inputchain->SetBranchAddress("TRKmeasurementsOnTrack_pixcl_sctcl_index",
0212                                  &TRKmeasurementsOnTrack_pixcl_sctcl_index);
0213   m_inputchain->SetBranchAddress("TRKoutliersOnTrack_pixcl_sctcl_index",
0214                                  &TRKoutliersOnTrack_pixcl_sctcl_index);
0215   m_inputchain->SetBranchAddress("TRKcharge", TRKcharge);
0216   m_inputchain->SetBranchAddress("TRKperigee_position", &TRKperigee_position);
0217   m_inputchain->SetBranchAddress("TRKperigee_momentum", &TRKperigee_momentum);
0218   m_inputchain->SetBranchAddress("TTCindex", TTCindex);
0219   m_inputchain->SetBranchAddress("TTCevent_index", TTCevent_index);
0220   m_inputchain->SetBranchAddress("TTCparticle_link", TTCparticle_link);
0221   m_inputchain->SetBranchAddress("TTCprobability", TTCprobability);
0222   m_inputchain->SetBranchAddress("nDTT", &nDTT);
0223   m_inputchain->SetBranchAddress("DTTindex", DTTindex);
0224   m_inputchain->SetBranchAddress("DTTsize", DTTsize);
0225   m_inputchain->SetBranchAddress("DTTtrajectory_eventindex",
0226                                  &DTTtrajectory_eventindex);
0227   m_inputchain->SetBranchAddress("DTTtrajectory_barcode",
0228                                  &DTTtrajectory_barcode);
0229   m_inputchain->SetBranchAddress("DTTstTruth_subDetType",
0230                                  &DTTstTruth_subDetType);
0231   m_inputchain->SetBranchAddress("DTTstTrack_subDetType",
0232                                  &DTTstTrack_subDetType);
0233   m_inputchain->SetBranchAddress("DTTstCommon_subDetType",
0234                                  &DTTstCommon_subDetType);
0235   */
0236
0237   for (const auto& file : m_cfg.inputfiles) {
0238     ACTS_DEBUG("Adding file '" << file << "' to tree " << m_cfg.treename);
0239     m_inputchain->Add(file.c_str());
0240   }
0241
0242   m_events = m_inputchain->GetEntries();
0243
0244   ACTS_DEBUG("End of constructor. In total available events=" << m_events);
0245 }  // constructor
0246
0247 SimParticleContainer RootAthenaDumpReader::readParticles() const {
0248   std::vector<SimParticle> particles;
0249   particles.reserve(nPartEVT);
0250
0251   for (auto ip = 0; ip < nPartEVT; ++ip) {
0252     if (m_cfg.onlyPassedParticles && !static_cast<bool>(Part_passed[ip])) {
0253       continue;
0254     }
0255
0256     SimBarcode dummyBarcode =
0257         SimBarcode()
0258             .withVertexPrimary(
0259                 static_cast<SimBarcode::PrimaryVertexId>(Part_event_number[ip]))
0260             .withVertexSecondary(static_cast<SimBarcode::SecondaryVertexId>(
0261                 Part_barcode[ip] < s_maxBarcodeForPrimary ? 0 : 1))
0262             .withParticle(
0263                 static_cast<SimBarcode::ParticleId>(Part_barcode[ip]));
0264     SimParticleState particle(dummyBarcode,
0265                               static_cast<Acts::PdgParticle>(Part_pdg_id[ip]));
0266
0267     Acts::Vector3 p = Acts::Vector3{Part_px[ip], Part_py[ip], Part_pz[ip]} *
0268                       Acts::UnitConstants::MeV;
0269     particle.setAbsoluteMomentum(p.norm());
0270
0271     particle.setDirection(p.normalized());
0272
0273     auto x = Acts::Vector4{Part_vx[ip], Part_vy[ip], Part_vz[ip], 0.0};
0274     particle.setPosition4(x);
0275
0276     particles.push_back(SimParticle(particle, particle));
0277   }
0278
0279   ACTS_DEBUG("Created " << particles.size() << " particles");
0280   auto before = particles.size();
0281
0282   auto particlesSet = particleVectorToSet(particles);
0283
0284   if (particlesSet.size() < before) {
0285     ACTS_WARNING("Particle IDs not unique for " << before - particles.size()
0286                                                 << " particles!");
0287   }
0288
0289   return particlesSet;
0290 }
0291
0292 std::tuple<ClusterContainer, MeasurementContainer, MeasurementParticlesMap,
0293            std::unordered_map<int, std::size_t>>
0294 RootAthenaDumpReader::readMeasurements(
0295     SimParticleContainer& particles, const Acts::GeometryContext& gctx) const {
0296   ClusterContainer clusters;
0297   clusters.reserve(nCL);
0298
0299   MeasurementContainer measurements;
0300   measurements.reserve(nCL);
0301
0302   std::size_t nTotalTotZero = 0;
0303
0304   const auto prevParticlesSize = particles.size();
0305   MeasurementParticlesMap measPartMap;
0306
0307   // We cannot use im for the index since we might skip measurements
0308   std::unordered_map<int, std::size_t> imIdxMap;
0309   imIdxMap.reserve(nCL);
0310
0311   for (int im = 0; im < nCL; im++) {
0312     if (!(CLhardware->at(im) == "PIXEL" || CLhardware->at(im) == "STRIP")) {
0313       ACTS_ERROR("hardware is neither 'PIXEL' or 'STRIP', skip particle");
0314       continue;
0315     }
0316     ACTS_VERBOSE("Cluster " << im << ": " << CLhardware->at(im));
0317
0318     auto type = (CLhardware->at(im) == "PIXEL") ? ePixel : eStrip;
0319
0320     // Make cluster
0321     // TODO refactor Cluster class so it is not so tedious
0322     Cluster cluster;
0323     if (m_cfg.readCellData && m_haveCellData) {
0324       const auto& etas = CLetas->at(im);
0325       const auto& phis = CLetas->at(im);
0326       const auto& tots = CLtots->at(im);
0327
0328       const auto totalTot = std::accumulate(tots.begin(), tots.end(), 0);
0329
0330       const auto [minEta, maxEta] =
0331           std::minmax_element(etas.begin(), etas.end());
0332       const auto [minPhi, maxPhi] =
0333           std::minmax_element(phis.begin(), phis.end());
0334
0335       cluster.channels.reserve(etas.size());
0336
0337       cluster.sizeLoc0 = *maxEta - *minEta;
0338       cluster.sizeLoc1 = *maxPhi - *minPhi;
0339
0340       if (totalTot == 0.0) {
0341         ACTS_VERBOSE(
0342             "total time over threshold is 0, set all activations to 0");
0343         nTotalTotZero++;
0344       }
0345
0346       for (const auto& [eta, phi, tot] : Acts::zip(etas, phis, tots)) {
0347         // Make best out of what we have:
0348         // Weight the overall collected charge corresponding to the
0349         // time-over-threshold of each cell Use this as activation (does this
0350         // make sense?)
0351         auto activation =
0352             (totalTot != 0.0) ? CLcharge_count[im] * tot / totalTot : 0.0;
0353
0354         // This bases every cluster at zero, but shouldn't matter right now
0355         ActsFatras::Segmentizer::Bin2D bin;
0356         bin[0] = eta - *minEta;
0357         bin[1] = phi - *minPhi;
0358
0359         // Of course we have no Segment2D because this is not Fatras
0360         cluster.channels.emplace_back(bin, ActsFatras::Segmentizer::Segment2D{},
0361                                       activation);
0362       }
0363
0364       ACTS_VERBOSE("- shape: " << cluster.channels.size()
0365                                << "cells, dimensions: " << cluster.sizeLoc0
0366                                << ", " << cluster.sizeLoc1);
0367     }
0368
0369     cluster.globalPosition = {CLx[im], CLy[im], CLz[im]};
0370     cluster.localDirection = {CLloc_direction1[im], CLloc_direction2[im],
0371                               CLloc_direction3[im]};
0372     cluster.lengthDirection = {CLJan_loc_direction1[im],
0373                                CLJan_loc_direction2[im],
0374                                CLJan_loc_direction3[im]};
0375     cluster.localEta = CLloc_eta[im];
0376     cluster.localPhi = CLloc_phi[im];
0377     cluster.globalEta = CLglob_eta[im];
0378     cluster.globalPhi = CLglob_phi[im];
0379     cluster.etaAngle = CLeta_angle[im];
0380     cluster.phiAngle = CLphi_angle[im];
0381
0382     // Measurement creation
0383     Acts::GeometryIdentifier geoId;
0384     std::vector<double> localParams;
0385     if (m_cfg.geometryIdMap && m_cfg.trackingGeometry) {
0386       const auto& geoIdMap = m_cfg.geometryIdMap->left;
0387       if (geoIdMap.find(CLmoduleID[im]) == geoIdMap.end()) {
0388         ACTS_WARNING("Missing geo id for " << CLmoduleID[im] << ", skip hit");
0389         continue;
0390       }
0391
0392       geoId = m_cfg.geometryIdMap->left.at(CLmoduleID[im]);
0393
0394       auto surface = m_cfg.trackingGeometry->findSurface(geoId);
0395       if (surface == nullptr) {
0396         ACTS_WARNING("Did not find " << geoId
0397                                      << " in tracking geometry, skip hit");
0398         continue;
0399       }
0400
0401       bool inside =
0402           surface->isOnSurface(gctx, cluster.globalPosition, {},
0403                                Acts::BoundaryTolerance::AbsoluteEuclidean(
0404                                    m_cfg.absBoundaryTolerance),
0405                                std::numeric_limits<double>::max());
0406
0407       if (!inside) {
0408         const Acts::Vector3 v =
0409             surface->localToGlobalTransform(gctx).inverse() *
0410             cluster.globalPosition;
0411         ACTS_WARNING("Projected position is not in surface bounds for "
0412                      << surface->geometryId() << ", skip hit");
0413         ACTS_WARNING("Position in local coordinates: " << v.transpose());
0414         ACTS_WARNING("Surface details:\n" << surface->toStream(gctx));
0415         continue;
0416       }
0417
0418       auto loc = surface->globalToLocal(gctx, cluster.globalPosition, {},
0419                                         Acts::s_onSurfaceTolerance);
0420
0421       if (!loc.ok()) {
0422         const Acts::Vector3 v =
0423             surface->localToGlobalTransform(gctx).inverse() *
0424             cluster.globalPosition;
0425         ACTS_WARNING("Global-to-local fit failed on "
0426                      << geoId << " (z dist: " << v[2]
0427                      << ", projected on surface: " << std::boolalpha << inside
0428                      << ") , skip hit");
0429         continue;
0430       }
0431
0432       // TODO is this in strip coordinates or in polar coordinates for annulus
0433       // bounds?
0434       localParams = std::vector<double>(loc->begin(), loc->end());
0435     } else {
0436       geoId = Acts::GeometryIdentifier(CLmoduleID[im]);
0437       localParams = {CLloc_direction1[im], CLloc_direction2[im]};
0438     }
0439
0440     const auto& locCov = CLlocal_cov->at(im);
0441
0442     DigitizedParameters digiPars;
0443     if (type == ePixel) {
0444       digiPars.indices = {Acts::eBoundLoc0, Acts::eBoundLoc1};
0445       assert(locCov.size() == 4);
0446       digiPars.variances = {locCov[0], locCov[3]};
0447       digiPars.values = localParams;
0448     } else {
0449       assert(!locCov.empty());
0450       // Barrel-endcap index can be -2/2 for endcap or 0 for barrel
0451       // We need to choose the coordinate of local measurement depending on that
0452       const static std::array boundLoc = {Acts::eBoundLoc0, Acts::eBoundLoc1};
0453       auto i = CLbarrel_endcap[im] == 0 ? 0 : 1;
0454       digiPars.variances = {locCov[i]};
0455       digiPars.values = {localParams[i]};
0456       digiPars.indices = {boundLoc[i]};
0457     }
0458
0459     std::size_t measIndex = measurements.size();
0460     ACTS_VERBOSE("Add measurement with index " << measIndex);
0461     imIdxMap.emplace(im, measIndex);
0462     createMeasurement(measurements, geoId, digiPars);
0463     cluster.geometryId = geoId;
0464     clusters.push_back(cluster);
0465
0466     if (!m_cfg.noTruth) {
0467       // Create measurement particles map and particles container
0468       for (const auto& [subevt, barcode] :
0469            Acts::zip(CLparticleLink_eventIndex->at(im),
0470                      CLparticleLink_barcode->at(im))) {
0471         SimBarcode dummyBarcode =
0472             SimBarcode()
0473                 .withVertexPrimary(
0474                     static_cast<SimBarcode::PrimaryVertexId>(subevt))
0475                 .withVertexSecondary(static_cast<SimBarcode::SecondaryVertexId>(
0476                     barcode < s_maxBarcodeForPrimary ? 0 : 1))
0477                 .withParticle(static_cast<SimBarcode::ParticleId>(barcode));
0478         // If we don't find the particle, create one with default values
0479         if (particles.find(dummyBarcode) == particles.end()) {
0480           ACTS_VERBOSE("Particle with subevt "
0481                        << subevt << ", barcode " << barcode
0482                        << "not found, create dummy one");
0483           particles.emplace(dummyBarcode, Acts::PdgParticle::eInvalid);
0484         }
0485         measPartMap.insert(
0486             std::pair<Index, SimBarcode>{measIndex, dummyBarcode});
0487       }
0488     }
0489   }
0490
0491   if (measurements.size() < static_cast<std::size_t>(nCL)) {
0492     ACTS_WARNING("Could not convert " << nCL - measurements.size() << " / "
0493                                       << nCL << " measurements");
0494   }
0495
0496   if (particles.size() - prevParticlesSize > 0) {
0497     ACTS_DEBUG("Created " << particles.size() - prevParticlesSize
0498                           << " dummy particles");
0499   }
0500
0501   if (nTotalTotZero > 0) {
0502     ACTS_DEBUG(nTotalTotZero << " / " << nCL
0503                              << " clusters have zero time-over-threshold");
0504   }
0505
0506   return {std::move(clusters), std::move(measurements), std::move(measPartMap),
0507           std::move(imIdxMap)};
0508 }
0509
0510 std::tuple<SpacePointContainer, SpacePointContainer, SpacePointContainer>
0511 RootAthenaDumpReader::readSpacePoints(
0512     const std::optional<std::unordered_map<int, std::size_t>>& imIdxMap) const {
0513   SpacePointContainer pixelSpacePoints(
0514       SpacePointColumns::SourceLinks | SpacePointColumns::X |
0515       SpacePointColumns::Y | SpacePointColumns::Z |
0516       SpacePointColumns::VarianceZ | SpacePointColumns::VarianceR);
0517   pixelSpacePoints.reserve(nSP);
0518
0519   SpacePointContainer stripSpacePoints(
0520       SpacePointColumns::SourceLinks | SpacePointColumns::X |
0521       SpacePointColumns::Y | SpacePointColumns::Z |
0522       SpacePointColumns::VarianceZ | SpacePointColumns::VarianceR |
0523       SpacePointColumns::StripCalibrationDetails);
0524   stripSpacePoints.reserve(nSP);
0525
0526   SpacePointContainer spacePoints(
0527       SpacePointColumns::SourceLinks | SpacePointColumns::X |
0528       SpacePointColumns::Y | SpacePointColumns::Z |
0529       SpacePointColumns::VarianceZ | SpacePointColumns::VarianceR |
0530       SpacePointColumns::StripCalibrationDetails);
0531   spacePoints.reserve(nSP);
0532
0533   // Loop on space points
0534   std::size_t skippedSpacePoints = 0;
0535   for (int isp = 0; isp < nSP; isp++) {
0536     auto isPhiOverlap = (SPisOverlap[isp] == 2) || (SPisOverlap[isp] == 3);
0537     auto isEtaOverlap = (SPisOverlap[isp] == 1) || (SPisOverlap[isp] == 3);
0538     if (m_cfg.skipOverlapSPsPhi && isPhiOverlap) {
0539       ++skippedSpacePoints;
0540       continue;
0541     }
0542     if (m_cfg.skipOverlapSPsEta && isEtaOverlap) {
0543       ++skippedSpacePoints;
0544       continue;
0545     }
0546
0547     const Acts::Vector3 globalPos{SPx[isp], SPy[isp], SPz[isp]};
0548     const double spCovr = SPcovr[isp];
0549     const double spCovz = SPcovz[isp];
0550
0551     // PIX=1  STRIP = 2
0552     auto type = SPCL2_index[isp] == -1 ? ePixel : eStrip;
0553
0554     ACTS_VERBOSE("SP:: " << type << " [" << globalPos.transpose() << "] "
0555                          << spCovr << " " << spCovz);
0556
0557     boost::container::static_vector<Acts::SourceLink, 2> sLinks;
0558
0559     const auto cl1Index = SPCL1_index[isp];
0560     assert(cl1Index >= 0 && cl1Index < nCL);
0561
0562     auto getGeoId =
0563         [&](auto athenaId) -> std::optional<Acts::GeometryIdentifier> {
0564       if (m_cfg.geometryIdMap == nullptr) {
0565         return Acts::GeometryIdentifier{athenaId};
0566       }
0567       if (m_cfg.geometryIdMap->left.find(athenaId) ==
0568           m_cfg.geometryIdMap->left.end()) {
0569         return std::nullopt;
0570       }
0571       return m_cfg.geometryIdMap->left.at(athenaId);
0572     };
0573
0574     auto cl1GeoId = getGeoId(CLmoduleID[cl1Index]);
0575     if (!cl1GeoId) {
0576       ACTS_WARNING("Could not find geoId for space point cluster 1");
0577       continue;
0578     }
0579
0580     if (imIdxMap && !imIdxMap->contains(cl1Index)) {
0581       ACTS_WARNING("Measurement 1 for space point " << isp << " not created");
0582       continue;
0583     }
0584
0585     IndexSourceLink first(*cl1GeoId,
0586                           imIdxMap ? imIdxMap->at(cl1Index) : cl1Index);
0587     sLinks.emplace_back(first);
0588
0589     // First create pixel space point here, later maybe overwrite with strip
0590     // space point
0591     auto sp = spacePoints.createSpacePoint();
0592     sp.x() = globalPos.x();
0593     sp.y() = globalPos.y();
0594     sp.z() = globalPos.z();
0595     sp.varianceR() = spCovr;
0596     sp.varianceZ() = spCovz;
0597
0598     if (type == ePixel) {
0599       auto pixelSp = pixelSpacePoints.createSpacePoint();
0600       pixelSp.assignSourceLinks(sLinks);
0601       pixelSp.x() = globalPos.x();
0602       pixelSp.y() = globalPos.y();
0603       pixelSp.z() = globalPos.z();
0604       pixelSp.varianceR() = spCovr;
0605       pixelSp.varianceZ() = spCovz;
0606     } else {
0607       const auto cl2Index = SPCL2_index[isp];
0608       assert(cl2Index >= 0 && cl2Index < nCL);
0609
0610       auto cl2GeoId = getGeoId(CLmoduleID[cl2Index]);
0611       if (!cl2GeoId) {
0612         ACTS_WARNING("Could not find geoId for space point cluster 2");
0613         continue;
0614       }
0615
0616       if (imIdxMap && !imIdxMap->contains(cl2Index)) {
0617         ACTS_WARNING("Measurement 2 for space point " << isp << " not created");
0618         continue;
0619       }
0620
0621       IndexSourceLink second(*cl2GeoId,
0622                              imIdxMap ? imIdxMap->at(cl2Index) : cl2Index);
0623       sLinks.emplace_back(second);
0624
0625       Acts::Vector3 outerStripCenter = Acts::Vector3::Zero();
0626       Acts::Vector3 stripSeparation = Acts::Vector3::Zero();
0627       Acts::Vector3 outerStripHalfVector = Acts::Vector3::Zero();
0628       Acts::Vector3 innerStripHalfVector = Acts::Vector3::Zero();
0629
0630       if (m_haveStripFeatures) {
0631         outerStripCenter = {SPtopStripCenterPosition->at(isp).at(0),
0632                             SPtopStripCenterPosition->at(isp).at(1),
0633                             SPtopStripCenterPosition->at(isp).at(2)};
0634         stripSeparation = {SPstripCenterDistance->at(isp).at(0),
0635                            SPstripCenterDistance->at(isp).at(1),
0636                            SPstripCenterDistance->at(isp).at(2)};
0637         outerStripHalfVector = {SPtopStripDirection->at(isp).at(0),
0638                                 SPtopStripDirection->at(isp).at(1),
0639                                 SPtopStripDirection->at(isp).at(2)};
0640         outerStripHalfVector *= SPhl_topstrip[isp];
0641         innerStripHalfVector = {SPbottomStripDirection->at(isp).at(0),
0642                                 SPbottomStripDirection->at(isp).at(1),
0643                                 SPbottomStripDirection->at(isp).at(2)};
0644         innerStripHalfVector *= SPhl_botstrip[isp];
0645       }
0646
0647       auto stripSp = stripSpacePoints.createSpacePoint();
0648       stripSp.assignSourceLinks(sLinks);
0649       stripSp.x() = globalPos.x();
0650       stripSp.y() = globalPos.y();
0651       stripSp.z() = globalPos.z();
0652       stripSp.varianceR() = spCovr;
0653       stripSp.varianceZ() = spCovz;
0654
0655       Eigen::Map<Eigen::Vector3f>(
0656           sp.outerStripCalibrationDetails().outerCenter.data()) =
0657           outerStripCenter.cast<float>();
0658       Eigen::Map<Eigen::Vector3f>(
0659           sp.outerStripCalibrationDetails().innerToOuterSeparation.data()) =
0660           stripSeparation.cast<float>();
0661       Eigen::Map<Eigen::Vector3f>(
0662           sp.outerStripCalibrationDetails().outerHalfVector.data()) =
0663           outerStripHalfVector.cast<float>();
0664       Eigen::Map<Eigen::Vector3f>(
0665           sp.outerStripCalibrationDetails().innerHalfVector.data()) =
0666           innerStripHalfVector.cast<float>();
0667
0668       sp.outerStripCalibrationDetails() =
0669           stripSp.outerStripCalibrationDetails();
0670     }
0671
0672     sp.assignSourceLinks(sLinks);
0673   }
0674
0675   if (m_cfg.skipOverlapSPsEta || m_cfg.skipOverlapSPsPhi) {
0676     ACTS_DEBUG("Skipped " << skippedSpacePoints
0677                           << " because of eta/phi overlaps");
0678   }
0679   if (spacePoints.size() <
0680       (static_cast<std::size_t>(nSP) - skippedSpacePoints)) {
0681     ACTS_WARNING("Could not convert " << nSP - spacePoints.size() << " of "
0682                                       << nSP << " space points");
0683   }
0684
0685   ACTS_DEBUG("Created " << spacePoints.size() << " overall space points");
0686   ACTS_DEBUG("Created " << pixelSpacePoints.size() << " "
0687                         << " pixel space points");
0688   ACTS_DEBUG("Created " << stripSpacePoints.size() << " "
0689                         << " strip space points");
0690
0691   return {std::move(spacePoints), std::move(pixelSpacePoints),
0692           std::move(stripSpacePoints)};
0693 }
0694
0695 std::pair<SimParticleContainer, MeasurementParticlesMap>
0696 RootAthenaDumpReader::reprocessParticles(
0697     const SimParticleContainer& particles,
0698     const MeasurementParticlesMap& measPartMap) const {
0699   std::vector<SimParticle> newParticles;
0700   newParticles.reserve(particles.size());
0701   MeasurementParticlesMap newMeasPartMap;
0702   newMeasPartMap.reserve(measPartMap.size());
0703
0704   const auto partMeasMap = invertIndexMultimap(measPartMap);
0705
0706   std::uint16_t primaryCount = 0;
0707   std::uint16_t secondaryCount = 0;
0708
0709   for (const auto& particle : particles) {
0710     const auto [begin, end] = partMeasMap.equal_range(particle.particleId());
0711
0712     if (begin == end) {
0713       ACTS_VERBOSE("Particle " << particle.particleId()
0714                                << " has no measurements");
0715       continue;
0716     }
0717
0718     auto primary = particle.particleId().vertexSecondary() == 0;
0719
0720     // vertex primary shouldn't be zero for a valid particle
0721     SimBarcode fatrasBarcode = SimBarcode().withVertexPrimary(1);
0722     if (primary) {
0723       fatrasBarcode =
0724           fatrasBarcode.withVertexSecondary(0).withParticle(primaryCount);
0725       assert(primaryCount < std::numeric_limits<std::uint16_t>::max());
0726       primaryCount++;
0727     } else {
0728       fatrasBarcode =
0729           fatrasBarcode.withVertexSecondary(1).withParticle(secondaryCount);
0730       assert(primaryCount < std::numeric_limits<std::uint16_t>::max());
0731       secondaryCount++;
0732     }
0733
0734     auto newParticle = particle.withParticleId(fatrasBarcode);
0735     newParticle.finalState().setNumberOfHits(std::distance(begin, end));
0736     newParticles.push_back(newParticle);
0737
0738     for (auto it = begin; it != end; ++it) {
0739       newMeasPartMap.insert(
0740           std::pair<Index, SimBarcode>{it->second, fatrasBarcode});
0741     }
0742   }
0743
0744   ACTS_DEBUG("After reprocessing particles " << newParticles.size() << " of "
0745                                              << particles.size() << " remain");
0746   return {particleVectorToSet(newParticles), std::move(newMeasPartMap)};
0747 }
0748
0749 ProcessCode RootAthenaDumpReader::read(const AlgorithmContext& ctx) {
0750   ACTS_DEBUG("Reading event " << ctx.eventNumber);
0751   auto entry = ctx.eventNumber;
0752   if (entry >= m_events) {
0753     ACTS_ERROR("event out of bounds");
0754     return ProcessCode::ABORT;
0755   }
0756
0757   std::lock_guard<std::mutex> lock(m_read_mutex);
0758
0759   m_inputchain->GetEntry(entry);
0760
0761   std::optional<std::unordered_map<int, std::size_t>> optImIdxMap;
0762
0763   if (!m_cfg.onlySpacePoints) {
0764     SimParticleContainer candidateParticles;
0765
0766     if (!m_cfg.noTruth) {
0767       candidateParticles = readParticles();
0768     }
0769
0770     auto [clusters, measurements, candidateMeasPartMap, imIdxMap] =
0771         readMeasurements(candidateParticles, ctx.simGeoContext);
0772     optImIdxMap.emplace(std::move(imIdxMap));
0773
0774     m_outputClusters(ctx, std::move(clusters));
0775     const auto& storedMeasurements =
0776         m_outputMeasurements(ctx, std::move(measurements));
0777     std::vector<MeasurementContainer::Index> allIndices(
0778         storedMeasurements.size());
0779     std::iota(allIndices.begin(), allIndices.end(),
0780               MeasurementContainer::Index{0});
0781     m_outputMeasurementSubset(
0782         ctx, MeasurementSubset(storedMeasurements, std::move(allIndices)));
0783
0784     if (!m_cfg.noTruth) {
0785       auto [particles, measPartMap] =
0786           reprocessParticles(candidateParticles, candidateMeasPartMap);
0787
0788       m_outputParticles(ctx, std::move(particles));
0789       m_outputParticleMeasMap(ctx, invertIndexMultimap(measPartMap));
0790       m_outputMeasParticleMap(ctx, std::move(measPartMap));
0791     }
0792   }
0793
0794   auto [spacePoints, pixelSpacePoints, stripSpacePoints] =
0795       readSpacePoints(optImIdxMap);
0796
0797   m_outputPixelSpacePoints(ctx, std::move(pixelSpacePoints));
0798   m_outputStripSpacePoints(ctx, std::move(stripSpacePoints));
0799   m_outputSpacePoints(ctx, std::move(spacePoints));
0800
0801   return ProcessCode::SUCCESS;
0802 }
0803 }  // namespace ActsExamples