Csv/src/CsvSeedWriter.cpp

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