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