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