Csv/src/CsvSeedWriter.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/Csv/CsvSeedWriter.hpp"
0010
0011 #include "ActsExamples/Utilities/EventDataTransforms.hpp"
0012 #include "ActsExamples/Utilities/Paths.hpp"
0013 #include "ActsExamples/Utilities/Range.hpp"
0014 #include "ActsExamples/Validation/TrackClassification.hpp"
0015
0016 #include <fstream>
0017 #include <ios>
0018 #include <iostream>
0019 #include <numbers>
0020 #include <stdexcept>
0021 #include <string>
0022 #include <unordered_map>
0023
0024 using Acts::VectorHelpers::phi;
0025 using Acts::VectorHelpers::theta;
0026
0027 namespace ActsExamples {
0028
0029 namespace {
0030
0031 /// @brief Struct for brief seed summary info
0032 ///
0033 struct SeedInfo {
0034   std::size_t seedID = 0;
0035   SimBarcode particleId;
0036   float seedPt = -1;
0037   float seedPhi = 0;
0038   float seedEta = 0;
0039   float vertexZ = 0;
0040   float quality = -1;
0041   boost::container::small_vector<Acts::Vector3, 3> globalPosition;
0042   float truthDistance = -1;
0043   std::string seedType = "unknown";
0044   ProtoTrack measurementsID;
0045 };
0046
0047 }  // namespace
0048
0049 CsvSeedWriter::CsvSeedWriter(const Config& config, Acts::Logging::Level level)
0050     : WriterT<TrackParametersContainer>(config.inputTrackParameters,
0051                                         "CsvSeedWriter", level),
0052       m_cfg(config) {
0053   if (m_cfg.inputSeeds.empty()) {
0054     throw std::invalid_argument("Missing space points input collection");
0055   }
0056   if (m_cfg.inputSimHits.empty()) {
0057     throw std::invalid_argument("Missing simulated hits input collection");
0058   }
0059   if (m_cfg.inputMeasurementParticlesMap.empty()) {
0060     throw std::invalid_argument("Missing hit-particles map input collection");
0061   }
0062   if (m_cfg.inputMeasurementSimHitsMap.empty()) {
0063     throw std::invalid_argument(
0064         "Missing hit-simulated-hits map input collection");
0065   }
0066   if (m_cfg.fileName.empty()) {
0067     throw std::invalid_argument("Missing output filename");
0068   }
0069   if (m_cfg.outputDir.empty()) {
0070     throw std::invalid_argument("Missing output directory");
0071   }
0072
0073   m_inputSeeds.initialize(m_cfg.inputSeeds);
0074   m_inputSimHits.initialize(m_cfg.inputSimHits);
0075   m_inputMeasurementParticlesMap.initialize(m_cfg.inputMeasurementParticlesMap);
0076   m_inputMeasurementSimHitsMap.initialize(m_cfg.inputMeasurementSimHitsMap);
0077 }
0078
0079 ProcessCode CsvSeedWriter::writeT(const AlgorithmContext& ctx,
0080                                   const TrackParametersContainer& trackParams) {
0081   // Read additional input collections
0082   const auto& seeds = m_inputSeeds(ctx);
0083   const auto& simHits = m_inputSimHits(ctx);
0084   const auto& hitParticlesMap = m_inputMeasurementParticlesMap(ctx);
0085   const auto& hitSimHitsMap = m_inputMeasurementSimHitsMap(ctx);
0086
0087   std::string path =
0088       perEventFilepath(m_cfg.outputDir, m_cfg.fileName, ctx.eventNumber);
0089
0090   std::ofstream mos(path, std::ofstream::out | std::ofstream::trunc);
0091   if (!mos) {
0092     throw std::ios_base::failure("Could not open '" + path + "' to write");
0093   }
0094
0095   std::unordered_map<std::size_t, SeedInfo> infoMap;
0096   std::unordered_map<SimBarcode, std::pair<std::size_t, float>> goodSeed;
0097
0098   // Loop over the estimated track parameters
0099   for (std::size_t iparams = 0; iparams < trackParams.size(); ++iparams) {
0100     // The estimated bound parameters vector
0101     const auto params = trackParams[iparams].parameters();
0102
0103     float seedPhi = params[Acts::eBoundPhi];
0104     float seedEta = std::atanh(std::cos(params[Acts::eBoundTheta]));
0105
0106     // Get the proto track from which the track parameters are estimated
0107     const auto& seed = seeds[iparams];
0108     const auto& ptrack = seedToProtoTrack(seed);
0109
0110     std::vector<ParticleHitCount> particleHitCounts;
0111     identifyContributingParticles(hitParticlesMap, ptrack, particleHitCounts);
0112     bool truthMatched = false;
0113     float truthDistance = -1;
0114     auto majorityParticleId = particleHitCounts.front().particleId;
0115     // Seed are considered truth matched if they have only one contributing
0116     // particle
0117     if (particleHitCounts.size() == 1) {
0118       truthMatched = true;
0119       // Get the index of the first space point
0120       const auto& hitIdx = ptrack.front();
0121       // Get the sim hits via the measurement to sim hits map
0122       auto indices = makeRange(hitSimHitsMap.equal_range(hitIdx));
0123       // Get the truth particle direction from the sim hits
0124       Acts::Vector3 truthUnitDir = {0, 0, 0};
0125       for (auto [_, simHitIdx] : indices) {
0126         const auto& simHit = *simHits.nth(simHitIdx);
0127         if (simHit.particleId() == majorityParticleId) {
0128           truthUnitDir = simHit.direction();
0129         }
0130       }
0131       // Compute the distance between the truth and estimated directions
0132       float truthPhi = phi(truthUnitDir);
0133       float truthEta = std::atanh(std::cos(theta(truthUnitDir)));
0134       float dEta = std::abs(truthEta - seedEta);
0135       float dPhi =
0136           std::abs(truthPhi - seedPhi) < std::numbers::pi_v<float>
0137               ? std::abs(truthPhi - seedPhi)
0138               : std::abs(truthPhi - seedPhi) - std::numbers::pi_v<float>;
0139       truthDistance = std::sqrt(dPhi * dPhi + dEta * dEta);
0140       // If the seed is truth matched, check if it is the closest one for the
0141       // contributing particle
0142       if (goodSeed.contains(majorityParticleId)) {
0143         if (goodSeed[majorityParticleId].second > truthDistance) {
0144           goodSeed[majorityParticleId] = std::make_pair(iparams, truthDistance);
0145         }
0146       } else {
0147         goodSeed[majorityParticleId] = std::make_pair(iparams, truthDistance);
0148       }
0149     }
0150     // Store the global position of the space points
0151     boost::container::small_vector<Acts::Vector3, 3> globalPosition;
0152     for (auto sp : seed.spacePoints()) {
0153       Acts::Vector3 pos(sp.x(), sp.y(), sp.z());
0154       globalPosition.push_back(pos);
0155     }
0156
0157     // track info
0158     SeedInfo toAdd;
0159     toAdd.seedID = iparams;
0160     toAdd.particleId = majorityParticleId;
0161     toAdd.seedPt = std::abs(1.0 / params[Acts::eBoundQOverP]) *
0162                    std::sin(params[Acts::eBoundTheta]);
0163     toAdd.seedPhi = seedPhi;
0164     toAdd.seedEta = seedEta;
0165     toAdd.vertexZ = seed.vertexZ();
0166     toAdd.quality = seed.quality();
0167     toAdd.globalPosition = globalPosition;
0168     toAdd.truthDistance = truthDistance;
0169     toAdd.seedType = truthMatched ? "duplicate" : "fake";
0170     toAdd.measurementsID = ptrack;
0171
0172     infoMap[toAdd.seedID] = toAdd;
0173   }
0174
0175   mos << "seed_id,particleId," << "pT,eta,phi," << "bX,bY,bZ," << "mX,mY,mZ,"
0176       << "tX,tY,tZ," << "good/duplicate/fake," << "vertexZ,quality,"
0177       << "Hits_ID" << '\n';
0178
0179   for (auto& [id, info] : infoMap) {
0180     if (goodSeed[info.particleId].first == id) {
0181       info.seedType = "good";
0182     }
0183     // write the track info
0184     mos << info.seedID << ",";
0185     mos << info.particleId << ",";
0186     mos << info.seedPt << ",";
0187     mos << info.seedEta << ",";
0188     mos << info.seedPhi << ",";
0189     for (auto& point : info.globalPosition) {
0190       mos << point.x() << ",";
0191       mos << point.y() << ",";
0192       mos << point.z() << ",";
0193     }
0194     mos << info.seedType << ",";
0195     mos << info.vertexZ << ",";
0196     mos << info.quality << ",";
0197     mos << "\"[";
0198     for (auto& ID : info.measurementsID) {
0199       mos << ID << ",";
0200     }
0201     mos << "]\"";
0202     mos << '\n';
0203   }
0204
0205   return ProcessCode::SUCCESS;
0206 }
0207
0208 }  // namespace ActsExamples