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