diff --git a/Examples/Framework/include/ActsExamples/EventData/Cluster.hpp b/Examples/Framework/include/ActsExamples/EventData/Cluster.hpp index d0e7bed55d6..3eb850b1214 100644 --- a/Examples/Framework/include/ActsExamples/EventData/Cluster.hpp +++ b/Examples/Framework/include/ActsExamples/EventData/Cluster.hpp @@ -25,6 +25,14 @@ struct Cluster { // TODO make this be provided by Fatras? Acts::Vector3 globalPosition = Acts::Vector3::Zero(); + Acts::Vector3 localDirection = Acts::Vector3::Zero(); + Acts::Vector3 lengthDirection = Acts::Vector3::Zero(); + float localEta = 0.f; + float localPhi = 0.f; + float globalEta = 0.f; + float globalPhi = 0.f; + float etaAngle = 0.f; + float phiAngle = 0.f; double sumActivations() const { return std::accumulate( diff --git a/Examples/Framework/include/ActsExamples/EventData/GeometryContainers.hpp b/Examples/Framework/include/ActsExamples/EventData/GeometryContainers.hpp index e6cec85a47f..0724ed8c4ae 100644 --- a/Examples/Framework/include/ActsExamples/EventData/GeometryContainers.hpp +++ b/Examples/Framework/include/ActsExamples/EventData/GeometryContainers.hpp @@ -20,6 +20,7 @@ #include #include +#include #include #include @@ -229,4 +230,9 @@ struct GeometryIdMultisetAccessor { const Container* container = nullptr; }; +/// A map that allows mapping back and forth between ACTS and Athena Geometry +/// Ids +using GeometryIdMapActsAthena = + boost::bimap; + } // namespace ActsExamples diff --git a/Examples/Io/Root/include/ActsExamples/Io/Root/RootAthenaDumpReader.hpp b/Examples/Io/Root/include/ActsExamples/Io/Root/RootAthenaDumpReader.hpp index 33919b4bae1..e917e07038d 100644 --- a/Examples/Io/Root/include/ActsExamples/Io/Root/RootAthenaDumpReader.hpp +++ b/Examples/Io/Root/include/ActsExamples/Io/Root/RootAthenaDumpReader.hpp @@ -9,6 +9,7 @@ #pragma once #include "Acts/Utilities/Logger.hpp" +#include "ActsExamples/EventData/GeometryContainers.hpp" #include "ActsExamples/EventData/Measurement.hpp" #include "ActsExamples/EventData/SimSpacePoint.hpp" #include "ActsExamples/Framework/DataHandle.hpp" @@ -62,13 +63,33 @@ class RootAthenaDumpReader : public IReader { // name of the track parameters (fitted by athena?) std::string outputTrackParameters = "athena_track_parameters"; + /// Only extract spacepoints + bool onlySpacepoints = false; + /// Only extract particles that passed the tracking requirements, for /// details see: /// https://gitlab.cern.ch/atlas/athena/-/blob/main/InnerDetector/InDetGNNTracking/src/DumpObjects.cxx?ref_type=heads#L1363 bool onlyPassedParticles = false; + /// Skip spacepoints with phi overlap bool skipOverlapSPsPhi = false; + + /// Skip spacepoints with eta overlap bool skipOverlapSPsEta = false; + + /// A map that provides a mapping between ACTS and Athena surface + /// identifiers + std::shared_ptr geometryIdMap = + nullptr; + + /// Tracking Geometry that contains the surfaces where we project + /// the measurements on + std::shared_ptr trackingGeometry = nullptr; + + /// When projecting measurements on ACTS surfaces, which euclidean boundary + /// tolerance should be allowed. If a value above zero is needed, this + /// indicates that the ACTS surfaces do not 100% include the athena surfaces + double absBoundaryTolerance = 0.0; }; RootAthenaDumpReader(const RootAthenaDumpReader &) = delete; @@ -94,12 +115,42 @@ class RootAthenaDumpReader : public IReader { const Config &config() const { return m_cfg; } private: + /// Particles with barcodes larger then this value are considered to be + /// secondary particles + /// https://gitlab.cern.ch/atlas/athena/-/blob/main/InnerDetector/InDetGNNTracking/src/DumpObjects.h?ref_type=heads#L101 + constexpr static int s_maxBarcodeForPrimary = 200000; + /// Private access to the logging instance const Acts::Logger &logger() const { return *m_logger; } /// The config class Config m_cfg; + /// Helper method to read particles + SimParticleContainer readParticles() const; + + /// Helper method to read measurements + std::tuple, + std::unordered_map> + readMeasurements(SimParticleContainer &particles, + const Acts::GeometryContext &gctx) const; + + /// Helper method to read spacepoints + /// @param imIdxMap optional remapping of indices. Since the measurement + /// index must be continuous, we need to remap the measurements indices + /// if we skip measurements in the first place + std::tuple + readSpacepoints(const std::optional> + &imIdxMap) const; + + /// Helper method to reprocess particle ids + std::pair> + reprocessParticles( + const SimParticleContainer &particles, + const IndexMultimap &measPartMap) const; + /// Write handlers WriteDataHandle m_outputPixelSpacePoints{ this, "outputPixelSpacepoints"}; diff --git a/Examples/Io/Root/src/RootAthenaDumpReader.cpp b/Examples/Io/Root/src/RootAthenaDumpReader.cpp index a76bee1997a..3ea44e78016 100644 --- a/Examples/Io/Root/src/RootAthenaDumpReader.cpp +++ b/Examples/Io/Root/src/RootAthenaDumpReader.cpp @@ -8,6 +8,7 @@ #include "ActsExamples/Io/Root/RootAthenaDumpReader.hpp" +#include "Acts/Definitions/Units.hpp" #include "Acts/EventData/SourceLink.hpp" #include "Acts/Geometry/GeometryIdentifier.hpp" #include "Acts/Utilities/Zip.hpp" @@ -21,60 +22,47 @@ #include #include -class BarcodeConstructor { - /// Particles with barcodes larger then this value are considered to be - /// secondary particles - /// https://gitlab.cern.ch/atlas/athena/-/blob/main/InnerDetector/InDetGNNTracking/src/DumpObjects.h?ref_type=heads#L101 - constexpr static int s_maxBarcodeForPrimary = 200000; - - std::uint16_t m_primaryCount = 0; - std::uint16_t m_secondaryCount = 0; - std::unordered_map m_barcodeMap; - - static std::uint64_t concatInts(int a, int b) { - auto va = static_cast(a); - auto vb = static_cast(b); - std::uint64_t value = (static_cast(va) << 32) | vb; - return value; - } - - public: - ActsFatras::Barcode getBarcode(int barcode, int evtnumber) { - auto v = concatInts(barcode, evtnumber); - auto found = m_barcodeMap.find(v); - if (found != m_barcodeMap.end()) { - return found->second; - } +using namespace Acts::UnitLiterals; - auto primary = (barcode < s_maxBarcodeForPrimary); +namespace { +std::uint64_t concatInts(int a, int b) { + auto va = static_cast(a); + auto vb = static_cast(b); + std::uint64_t value = (static_cast(va) << 32) | vb; + return value; +} - ActsFatras::Barcode fBarcode; +std::pair splitInt(std::uint64_t v) { + return {static_cast((v & 0xFFFFFFFF00000000LL) >> 32), + static_cast(v & 0xFFFFFFFFLL)}; +} - // vertex primary shouldn't be zero for a valid particle - fBarcode.setVertexPrimary(1); - if (primary) { - fBarcode.setVertexSecondary(0); - fBarcode.setParticle(m_primaryCount); - assert(m_primaryCount < std::numeric_limits::max()); - m_primaryCount++; - } else { - fBarcode.setVertexSecondary(1); - fBarcode.setParticle(m_secondaryCount); - assert(m_primaryCount < std::numeric_limits::max()); - m_secondaryCount++; - } +/// In cases when there is built up a particle collection in an iterative way it +/// can be way faster to build up a vector and afterwards use a special +/// constructor to speed up the set creation. +inline auto particleVectorToSet(std::vector& particles) { + using namespace ActsExamples; + auto cmp = [](const auto& a, const auto& b) { + return a.particleId().value() == b.particleId().value(); + }; + + std::sort(particles.begin(), particles.end(), detail::CompareParticleId{}); + particles.erase(std::unique(particles.begin(), particles.end(), cmp), + particles.end()); + + return SimParticleContainer(boost::container::ordered_unique_range_t{}, + particles.begin(), particles.end()); +} - m_barcodeMap[v] = fBarcode; - return fBarcode; - } -}; +} // namespace enum SpacePointType { ePixel = 1, eStrip = 2 }; -ActsExamples::RootAthenaDumpReader::RootAthenaDumpReader( - const ActsExamples::RootAthenaDumpReader::Config& config, - Acts::Logging::Level level) - : ActsExamples::IReader(), +namespace ActsExamples { + +RootAthenaDumpReader::RootAthenaDumpReader( + const RootAthenaDumpReader::Config& config, Acts::Logging::Level level) + : IReader(), m_cfg(config), m_logger(Acts::getDefaultLogger(name(), level)) { if (m_cfg.inputfile.empty()) { @@ -89,10 +77,12 @@ ActsExamples::RootAthenaDumpReader::RootAthenaDumpReader( m_outputPixelSpacePoints.initialize(m_cfg.outputPixelSpacePoints); m_outputStripSpacePoints.initialize(m_cfg.outputStripSpacePoints); m_outputSpacePoints.initialize(m_cfg.outputSpacePoints); - m_outputClusters.initialize(m_cfg.outputClusters); - m_outputParticles.initialize(m_cfg.outputParticles); - m_outputMeasParticleMap.initialize(m_cfg.outputMeasurementParticlesMap); - m_outputMeasurements.initialize(m_cfg.outputMeasurements); + if (!m_cfg.onlySpacepoints) { + m_outputClusters.initialize(m_cfg.outputClusters); + m_outputParticles.initialize(m_cfg.outputParticles); + m_outputMeasParticleMap.initialize(m_cfg.outputMeasurementParticlesMap); + m_outputMeasurements.initialize(m_cfg.outputMeasurements); + } // Set the branches @@ -252,35 +242,19 @@ ActsExamples::RootAthenaDumpReader::RootAthenaDumpReader( m_events = m_inputchain->GetEntries(); ACTS_DEBUG("End of constructor. In total available events=" << m_events); - } // constructor -ActsExamples::ProcessCode ActsExamples::RootAthenaDumpReader::read( - const ActsExamples::AlgorithmContext& ctx) { - ACTS_DEBUG("Reading event " << ctx.eventNumber); - auto entry = ctx.eventNumber; - if (entry >= m_events) { - ACTS_ERROR("event out of bounds"); - return ProcessCode::ABORT; - } - - std::lock_guard lock(m_read_mutex); - - m_inputchain->GetEntry(entry); - - // Concat the two 32bit integers from athena to a Fatras barcode - - SimParticleContainer particles; - BarcodeConstructor barcodeConstructor; +SimParticleContainer RootAthenaDumpReader::readParticles() const { + std::vector particles; + particles.reserve(nPartEVT); for (auto ip = 0; ip < nPartEVT; ++ip) { if (m_cfg.onlyPassedParticles && !static_cast(Part_passed[ip])) { continue; } - auto barcode = - barcodeConstructor.getBarcode(Part_barcode[ip], Part_event_number[ip]); - SimParticle particle(barcode, + auto dummyBarcode = concatInts(Part_barcode[ip], Part_event_number[ip]); + SimParticle particle(dummyBarcode, static_cast(Part_pdg_id[ip])); Acts::Vector3 p = Acts::Vector3{Part_px[ip], Part_py[ip], Part_pz[ip]} * @@ -292,26 +266,53 @@ ActsExamples::ProcessCode ActsExamples::RootAthenaDumpReader::read( auto x = Acts::Vector4{Part_vx[ip], Part_vy[ip], Part_vz[ip], 0.0}; particle.setPosition4(x); - particles.insert(particle); + particles.push_back(particle); } + ACTS_DEBUG("Created " << particles.size() << " particles"); + auto before = particles.size(); + + auto particlesSet = particleVectorToSet(particles); + + if (particlesSet.size() < before) { + ACTS_WARNING("Particle IDs not unique for " << before - particles.size() + << " particles!"); + } + + return particlesSet; +} + +std::tuple, + std::unordered_map> +RootAthenaDumpReader::readMeasurements( + SimParticleContainer& particles, const Acts::GeometryContext& gctx) const { ClusterContainer clusters(nCL); + clusters.reserve(nCL); + MeasurementContainer measurements; measurements.reserve(nCL); + std::size_t nTotalTotZero = 0; + + const auto prevParticlesSize = particles.size(); IndexMultimap measPartMap; + // We cannot use im for the index since we might skip measurements + std::size_t idx = 0; + std::unordered_map imIdxMap; + for (int im = 0; im < nCL; im++) { if (!(CLhardware->at(im) == "PIXEL" || CLhardware->at(im) == "STRIP")) { - ACTS_ERROR("hardware is neither 'PIXEL' or 'STRIP'"); - return ActsExamples::ProcessCode::ABORT; + ACTS_ERROR("hardware is neither 'PIXEL' or 'STRIP', skip particle"); + continue; } ACTS_VERBOSE("Cluster " << im << ": " << CLhardware->at(im)); auto type = (CLhardware->at(im) == "PIXEL") ? ePixel : eStrip; // Make cluster - // TODO refactor ActsExamples::Cluster class so it is not so tedious + // TODO refactor Cluster class so it is not so tedious Cluster cluster; const auto& etas = CLetas->at(im); @@ -326,12 +327,18 @@ ActsExamples::ProcessCode ActsExamples::RootAthenaDumpReader::read( cluster.sizeLoc0 = *maxEta - *minEta; cluster.sizeLoc1 = *maxPhi - *minPhi; + if (totalTot == 0.0) { + ACTS_VERBOSE("total time over threshold is 0, set all activations to 0"); + nTotalTotZero++; + } + for (const auto& [eta, phi, tot] : Acts::zip(etas, phis, tots)) { // Make best out of what we have: // Weight the overall collected charge corresponding to the // time-over-threshold of each cell Use this as activation (does this make // sense?) - auto activation = CLcharge_count[im] * tot / totalTot; + auto activation = + (totalTot != 0.0) ? CLcharge_count[im] * tot / totalTot : 0.0; // This bases every cluster at zero, but shouldn't matter right now ActsFatras::Segmentizer::Bin2D bin; @@ -344,6 +351,17 @@ ActsExamples::ProcessCode ActsExamples::RootAthenaDumpReader::read( } cluster.globalPosition = {CLx[im], CLy[im], CLz[im]}; + cluster.localDirection = {CLloc_direction1[im], CLloc_direction2[im], + CLloc_direction3[im]}; + cluster.lengthDirection = {CLJan_loc_direction1[im], + CLJan_loc_direction2[im], + CLJan_loc_direction3[im]}; + cluster.localEta = CLloc_eta[im]; + cluster.localPhi = CLloc_phi[im]; + cluster.globalEta = CLglob_eta[im]; + cluster.globalPhi = CLglob_phi[im]; + cluster.etaAngle = CLeta_angle[im]; + cluster.phiAngle = CLphi_angle[im]; ACTS_VERBOSE("CL shape: " << cluster.channels.size() << "cells, dimensions: " << cluster.sizeLoc0 @@ -357,37 +375,128 @@ ActsExamples::ProcessCode ActsExamples::RootAthenaDumpReader::read( << CLloc_direction3[im]); const auto& locCov = CLlocal_cov->at(im); + std::optional sl; + std::vector localParams; + if (m_cfg.geometryIdMap && m_cfg.trackingGeometry) { + const auto& geoIdMap = m_cfg.geometryIdMap->left; + if (geoIdMap.find(CLmoduleID[im]) == geoIdMap.end()) { + ACTS_WARNING("Missing geo id for " << CLmoduleID[im] << ", skip hit"); + continue; + } + + auto geoId = m_cfg.geometryIdMap->left.at(CLmoduleID[im]); + sl = IndexSourceLink(geoId, idx); + + auto surface = m_cfg.trackingGeometry->findSurface(geoId); + if (surface == nullptr) { + ACTS_WARNING("Did not find " << geoId + << " in tracking geometry, skip hit"); + continue; + } + + bool inside = + surface->isOnSurface(gctx, cluster.globalPosition, {}, + Acts::BoundaryTolerance::AbsoluteEuclidean{ + m_cfg.absBoundaryTolerance}, + std::numeric_limits::max()); + + if (!inside) { + const Acts::Vector3 v = + surface->transform(gctx).inverse() * cluster.globalPosition; + ACTS_WARNING("Projected position is not in surface bounds for " + << surface->geometryId() << ", skip hit"); + ACTS_WARNING("Position in local coordinates: " << v.transpose()); + ACTS_WARNING("Surface details:\n" << surface->toStream(gctx)); + continue; + } + + const double tol = (type == ePixel) ? Acts::s_onSurfaceTolerance : 1.3_mm; + auto loc = surface->globalToLocal(gctx, cluster.globalPosition, {}, tol); + + if (!loc.ok()) { + const Acts::Vector3 v = + surface->transform(gctx).inverse() * cluster.globalPosition; + ACTS_WARNING("Global-to-local fit failed on " + << geoId << " (z dist: " << v[2] + << ", projected on surface: " << std::boolalpha << inside + << ") , skip hit"); + continue; + } + + // TODO is this in strip coordinates or in polar coordinates for annulus + // bounds? + localParams = std::vector(loc->begin(), loc->end()); + } else { + sl = IndexSourceLink(Acts::GeometryIdentifier(CLmoduleID[im]), idx); + localParams = {CLloc_direction1[im], CLloc_direction2[im]}; + } + DigitizedParameters digiPars; if (type == ePixel) { digiPars.indices = {Acts::eBoundLoc0, Acts::eBoundLoc1}; - digiPars.values = {CLloc_direction1[im], CLloc_direction2[im]}; assert(locCov.size() == 4); digiPars.variances = {locCov[0], locCov[3]}; + digiPars.values = localParams; } else { - digiPars.values = {CLloc_direction1[im]}; digiPars.indices = {Acts::eBoundLoc0}; assert(!locCov.empty()); digiPars.variances = {locCov[0]}; + digiPars.values = {localParams[0]}; } - IndexSourceLink sl(Acts::GeometryIdentifier{CLmoduleID[im]}, im); - - createMeasurement(measurements, digiPars, sl); + createMeasurement(measurements, digiPars, *sl); // Create measurement particles map and particles container - for (const auto& [subevt, bc] : Acts::zip(CLparticleLink_eventIndex->at(im), - CLparticleLink_barcode->at(im))) { - auto barcode = barcodeConstructor.getBarcode(bc, subevt); - measPartMap.insert(std::pair{im, barcode}); + for (const auto& [subevt, barcode] : + Acts::zip(CLparticleLink_eventIndex->at(im), + CLparticleLink_barcode->at(im))) { + auto dummyBarcode = concatInts(barcode, subevt); + // If we don't find the particle, create one with default values + if (particles.find(dummyBarcode) == particles.end()) { + ACTS_VERBOSE("Particle with subevt " << subevt << ", barcode " + << barcode + << "not found, create dummy one"); + particles.emplace(dummyBarcode, Acts::PdgParticle::eInvalid); + } + measPartMap.insert( + std::pair{idx, dummyBarcode}); } + + // Finally increment the measurement index + imIdxMap.emplace(im, idx); + ++idx; + } + + if (measurements.size() < static_cast(nCL)) { + ACTS_WARNING("Could not convert " << nCL - measurements.size() << " / " + << nCL << " measurements"); + } + + if (particles.size() - prevParticlesSize > 0) { + ACTS_DEBUG("Created " << particles.size() - prevParticlesSize + << " dummy particles"); + } + + if (nTotalTotZero > 0) { + ACTS_WARNING(nTotalTotZero << " / " << nCL + << " clusters have zero time-over-threshold"); } - // Prepare pixel space points + return {clusters, measurements, measPartMap, imIdxMap}; +} + +std::tuple +RootAthenaDumpReader::readSpacepoints( + const std::optional>& imIdxMap) const { SimSpacePointContainer pixelSpacePoints; + pixelSpacePoints.reserve(nSP); + + SimSpacePointContainer stripSpacePoints; + stripSpacePoints.reserve(nSP); - // Prepare space-point container - // They contain both pixel and SCT space points SimSpacePointContainer spacePoints; + spacePoints.reserve(nSP); // Loop on space points std::size_t skippedSpacePoints = 0; @@ -403,63 +512,211 @@ ActsExamples::ProcessCode ActsExamples::RootAthenaDumpReader::read( continue; } - Acts::Vector3 globalPos{SPx[isp], SPy[isp], SPz[isp]}; - double sp_covr = SPcovr[isp]; - double sp_covz = SPcovz[isp]; + const Acts::Vector3 globalPos{SPx[isp], SPy[isp], SPz[isp]}; + const double spCovr = SPcovr[isp]; + const double spCovz = SPcovz[isp]; // PIX=1 STRIP = 2 auto type = SPCL2_index[isp] == -1 ? ePixel : eStrip; ACTS_VERBOSE("SP:: " << type << " [" << globalPos.transpose() << "] " - << sp_covr << " " << sp_covz); + << spCovr << " " << spCovz); boost::container::static_vector sLinks; const auto cl1Index = SPCL1_index[isp]; assert(cl1Index >= 0 && cl1Index < nCL); - // NOTE This of course does not produce a valid Acts-stlye geometry id, but - // we can use it for the module map - IndexSourceLink first(Acts::GeometryIdentifier{CLmoduleID[cl1Index]}, - cl1Index); - sLinks.emplace_back(first); - - if (type == eStrip) { - const auto cl2Index = SPCL2_index[isp]; - assert(cl2Index >= 0 && cl2Index < nCL); + auto getGeoId = + [&](auto athenaId) -> std::optional { + if (m_cfg.geometryIdMap == nullptr) { + return Acts::GeometryIdentifier{athenaId}; + } + if (m_cfg.geometryIdMap->left.find(athenaId) == + m_cfg.geometryIdMap->left.end()) { + return std::nullopt; + } + return m_cfg.geometryIdMap->left.at(athenaId); + }; + + auto cl1GeoId = getGeoId(CLmoduleID[cl1Index]); + if (!cl1GeoId) { + ACTS_WARNING("Could not find geoId for spacepoint cluster 1"); + continue; + } - // NOTE This of course does not produce a valid Acts-stlye geometry id, - // but we can use it for the module map - IndexSourceLink second(Acts::GeometryIdentifier{CLmoduleID[cl2Index]}, - cl2Index); - sLinks.emplace_back(second); + if (imIdxMap && !imIdxMap->contains(cl1Index)) { + ACTS_WARNING("Measurement 1 for spacepoint " << isp << " not created"); + continue; } - SimSpacePoint sp(globalPos, std::nullopt, sp_covr, sp_covz, std::nullopt, + IndexSourceLink first(*cl1GeoId, + imIdxMap ? imIdxMap->at(cl1Index) : cl1Index); + sLinks.emplace_back(first); + + // First create pixel spacepoint here, later maybe overwrite with strip + // spacepoint + SimSpacePoint sp(globalPos, std::nullopt, spCovr, spCovz, std::nullopt, sLinks); if (type == ePixel) { pixelSpacePoints.push_back(sp); + } else { + const auto cl2Index = SPCL2_index[isp]; + assert(cl2Index >= 0 && cl2Index < nCL); + + auto cl2GeoId = getGeoId(CLmoduleID[cl1Index]); + if (!cl2GeoId) { + ACTS_WARNING("Could not find geoId for spacepoint cluster 2"); + continue; + } + + if (imIdxMap && !imIdxMap->contains(cl2Index)) { + ACTS_WARNING("Measurement 2 for spacepoint " << isp << " not created"); + continue; + } + + IndexSourceLink second(*cl2GeoId, + imIdxMap ? imIdxMap->at(cl2Index) : cl2Index); + sLinks.emplace_back(second); + + using Vector3f = Eigen::Matrix; + const Vector3f topStripDirection{SPtopStripDirection->at(isp).at(0), + SPtopStripDirection->at(isp).at(1), + SPtopStripDirection->at(isp).at(2)}; + const Vector3f bottomStripDirection{ + SPbottomStripDirection->at(isp).at(0), + SPbottomStripDirection->at(isp).at(1), + SPbottomStripDirection->at(isp).at(2)}; + const Vector3f stripCenterDistance{SPstripCenterDistance->at(isp).at(0), + SPstripCenterDistance->at(isp).at(1), + SPstripCenterDistance->at(isp).at(2)}; + const Vector3f topStripCenterPosition{ + SPtopStripCenterPosition->at(isp).at(0), + SPtopStripCenterPosition->at(isp).at(1), + SPtopStripCenterPosition->at(isp).at(2)}; + + sp = SimSpacePoint(globalPos, std::nullopt, spCovr, spCovz, std::nullopt, + sLinks, SPhl_topstrip[isp], SPhl_botstrip[isp], + topStripDirection.cast(), + bottomStripDirection.cast(), + stripCenterDistance.cast(), + topStripCenterPosition.cast()); + + stripSpacePoints.push_back(sp); } spacePoints.push_back(sp); } - ACTS_DEBUG("Created " << particles.size() << " particles"); if (m_cfg.skipOverlapSPsEta || m_cfg.skipOverlapSPsPhi) { ACTS_DEBUG("Skipped " << skippedSpacePoints << " because of eta/phi overlaps"); } + if (spacePoints.size() < static_cast(nSP)) { + ACTS_WARNING("Could not convert " << nSP - spacePoints.size() << " of " + << nSP << " spacepoints"); + } + ACTS_DEBUG("Created " << spacePoints.size() << " overall space points"); ACTS_DEBUG("Created " << pixelSpacePoints.size() << " " << " pixel space points"); + return {spacePoints, pixelSpacePoints, stripSpacePoints}; +} + +std::pair> +RootAthenaDumpReader::reprocessParticles( + const SimParticleContainer& particles, + const IndexMultimap& measPartMap) const { + std::vector newParticles; + newParticles.reserve(particles.size()); + IndexMultimap newMeasPartMap; + + const auto partMeasMap = invertIndexMultimap(measPartMap); + + std::uint16_t primaryCount = 0; + std::uint16_t secondaryCount = 0; + + for (const auto& particle : particles) { + const auto [begin, end] = partMeasMap.equal_range(particle.particleId()); + + if (begin == end) { + ACTS_VERBOSE("Particle " << particle.particleId().value() + << " has no measurements"); + continue; + } + + auto [athBarcode, athSubevent] = splitInt(particle.particleId().value()); + auto primary = (athBarcode < s_maxBarcodeForPrimary); + + ActsFatras::Barcode fatrasBarcode; + + // vertex primary shouldn't be zero for a valid particle + fatrasBarcode.setVertexPrimary(1); + if (primary) { + fatrasBarcode.setVertexSecondary(0); + fatrasBarcode.setParticle(primaryCount); + assert(primaryCount < std::numeric_limits::max()); + primaryCount++; + } else { + fatrasBarcode.setVertexSecondary(1); + fatrasBarcode.setParticle(secondaryCount); + assert(primaryCount < std::numeric_limits::max()); + secondaryCount++; + } + + auto newParticle = particle.withParticleId(fatrasBarcode); + newParticles.push_back(newParticle); + + for (auto it = begin; it != end; ++it) { + newMeasPartMap.insert( + std::pair{it->second, fatrasBarcode}); + } + } + + ACTS_DEBUG("After reprocessing particles " << newParticles.size() << " of " + << particles.size() << " remain"); + return {particleVectorToSet(newParticles), newMeasPartMap}; +} + +ProcessCode RootAthenaDumpReader::read(const AlgorithmContext& ctx) { + ACTS_DEBUG("Reading event " << ctx.eventNumber); + auto entry = ctx.eventNumber; + if (entry >= m_events) { + ACTS_ERROR("event out of bounds"); + return ProcessCode::ABORT; + } + + std::lock_guard lock(m_read_mutex); + + m_inputchain->GetEntry(entry); + + std::optional> optImIdxMap; + + if (!m_cfg.onlySpacepoints) { + auto candidateParticles = readParticles(); + + auto [clusters, measurements, candidateMeasPartMap, imIdxMap] = + readMeasurements(candidateParticles, ctx.geoContext); + optImIdxMap.emplace(std::move(imIdxMap)); + + auto [particles, measPartMap] = + reprocessParticles(candidateParticles, candidateMeasPartMap); + + m_outputClusters(ctx, std::move(clusters)); + m_outputParticles(ctx, std::move(particles)); + m_outputMeasParticleMap(ctx, std::move(measPartMap)); + m_outputMeasurements(ctx, std::move(measurements)); + } + + auto [spacePoints, pixelSpacePoints, stripSpacePoints] = + readSpacepoints(optImIdxMap); + m_outputPixelSpacePoints(ctx, std::move(pixelSpacePoints)); + m_outputStripSpacePoints(ctx, std::move(stripSpacePoints)); m_outputSpacePoints(ctx, std::move(spacePoints)); - m_outputClusters(ctx, std::move(clusters)); - m_outputParticles(ctx, std::move(particles)); - m_outputMeasParticleMap(ctx, std::move(measPartMap)); - m_outputMeasurements(ctx, std::move(measurements)); return ProcessCode::SUCCESS; } +} // namespace ActsExamples diff --git a/Examples/Python/src/Input.cpp b/Examples/Python/src/Input.cpp index f75af4f4271..e876efaa69c 100644 --- a/Examples/Python/src/Input.cpp +++ b/Examples/Python/src/Input.cpp @@ -93,7 +93,16 @@ void addInput(Context& ctx) { ACTS_PYTHON_DECLARE_READER( ActsExamples::RootAthenaDumpReader, mex, "RootAthenaDumpReader", treename, inputfile, outputMeasurements, outputPixelSpacePoints, - outputStripSpacePoints, outputSpacePoints, outputClusters); + outputStripSpacePoints, outputSpacePoints, outputClusters, + outputMeasurementParticlesMap, outputParticles, onlyPassedParticles, + skipOverlapSPsPhi, skipOverlapSPsEta, geometryIdMap, trackingGeometry, + absBoundaryTolerance); + +#ifdef WITH_GEOMODEL_PLUGIN + ACTS_PYTHON_DECLARE_READER(ActsExamples::RootAthenaDumpGeoIdCollector, mex, + "RootAthenaDumpGeoIdCollector", treename, + inputfile, trackingGeometry, geometryIdMap); +#endif ACTS_PYTHON_DECLARE_READER(ActsExamples::RootSimHitReader, mex, "RootSimHitReader", treeName, filePath,