diff --git a/Examples/Io/Root/include/ActsExamples/Io/Root/RootTrackParameterWriter.hpp b/Examples/Io/Root/include/ActsExamples/Io/Root/RootTrackParameterWriter.hpp index a54828983ad..625ef34ace0 100644 --- a/Examples/Io/Root/include/ActsExamples/Io/Root/RootTrackParameterWriter.hpp +++ b/Examples/Io/Root/include/ActsExamples/Io/Root/RootTrackParameterWriter.hpp @@ -1,6 +1,6 @@ // This file is part of the Acts project. // -// Copyright (C) 2017-2018 CERN for the benefit of the Acts project +// Copyright (C) 2017-2024 CERN for the benefit of the Acts project // // This Source Code Form is subject to the terms of the Mozilla Public // License, v. 2.0. If a copy of the MPL was not distributed with this @@ -23,9 +23,6 @@ class TFile; class TTree; -namespace ActsFatras { -class Barcode; -} // namespace ActsFatras namespace ActsExamples { struct AlgorithmContext; @@ -95,29 +92,68 @@ class RootTrackParameterWriter final : public TrackParameterWriter { ReadDataHandle m_inputMeasurementSimHitsMap{ this, "InputMeasurementSimHitsMap"}; - std::mutex m_writeMutex; ///< Mutex used to protect multi-threaded writes - TFile* m_outputFile{nullptr}; ///< The output file - TTree* m_outputTree{nullptr}; ///< The output tree - int m_eventNr{0}; ///< the event number of - - float m_loc0{NaNfloat}; ///< loc0 - float m_loc1{NaNfloat}; ///< loc1 - float m_phi{NaNfloat}; ///< phi - float m_theta{NaNfloat}; ///< theta - float m_qop{NaNfloat}; ///< q/p - float m_time{NaNfloat}; ///< time - float m_p{NaNfloat}; ///< p - float m_pt{NaNfloat}; ///< pt - float m_eta{NaNfloat}; ///< eta - - int m_t_charge{0}; ///< Truth particle charge - float m_t_loc0{NaNfloat}; ///< Truth parameter loc0 - float m_t_loc1{NaNfloat}; ///< Truth parameter loc1 - float m_t_phi{NaNfloat}; ///< Truth parameter phi - float m_t_theta{NaNfloat}; ///< Truth parameter theta - float m_t_qop{NaNfloat}; ///< Truth parameter qop - float m_t_time{NaNfloat}; ///< Truth parameter time - bool m_truthMatched = false; ///< Whether the seed is matched with truth + /// Mutex used to protect multi-threaded writes + std::mutex m_writeMutex; + TFile* m_outputFile{nullptr}; + TTree* m_outputTree{nullptr}; + + int m_eventNr{0}; + + int m_volumeId{0}; + int m_layerId{0}; + int m_surfaceId{0}; + + // Track parameters + float m_loc0{NaNfloat}; + float m_loc1{NaNfloat}; + float m_phi{NaNfloat}; + float m_theta{NaNfloat}; + float m_qop{NaNfloat}; + float m_time{NaNfloat}; + + float m_err_loc0{NaNfloat}; + float m_err_loc1{NaNfloat}; + float m_err_phi{NaNfloat}; + float m_err_theta{NaNfloat}; + float m_err_qop{NaNfloat}; + float m_err_time{NaNfloat}; + + int m_charge{0}; + float m_p{NaNfloat}; + float m_pt{NaNfloat}; + float m_eta{NaNfloat}; + + // Truth parameters + /// Whether the seed is matched with truth + bool m_t_matched{false}; + std::uint64_t m_t_particleId{0}; + unsigned int m_nMajorityHits{0}; + + float m_t_loc0{NaNfloat}; + float m_t_loc1{NaNfloat}; + float m_t_phi{NaNfloat}; + float m_t_theta{NaNfloat}; + float m_t_qop{NaNfloat}; + float m_t_time{NaNfloat}; + + int m_t_charge{0}; + float m_t_p{NaNfloat}; + float m_t_pt{NaNfloat}; + float m_t_eta{NaNfloat}; + + float m_res_loc0{NaNfloat}; + float m_res_loc1{NaNfloat}; + float m_res_phi{NaNfloat}; + float m_res_theta{NaNfloat}; + float m_res_qop{NaNfloat}; + float m_res_time{NaNfloat}; + + float m_pull_loc0{NaNfloat}; + float m_pull_loc1{NaNfloat}; + float m_pull_phi{NaNfloat}; + float m_pull_theta{NaNfloat}; + float m_pull_qop{NaNfloat}; + float m_pull_time{NaNfloat}; }; } // namespace ActsExamples diff --git a/Examples/Io/Root/src/RootTrackParameterWriter.cpp b/Examples/Io/Root/src/RootTrackParameterWriter.cpp index e7ea5f909af..57007ef64d8 100644 --- a/Examples/Io/Root/src/RootTrackParameterWriter.cpp +++ b/Examples/Io/Root/src/RootTrackParameterWriter.cpp @@ -13,6 +13,7 @@ #include "Acts/Utilities/MultiIndex.hpp" #include "ActsExamples/EventData/AverageSimHits.hpp" #include "ActsExamples/Framework/AlgorithmContext.hpp" +#include "ActsExamples/Framework/WriterT.hpp" #include "ActsExamples/Utilities/Range.hpp" #include "ActsExamples/Validation/TrackClassification.hpp" #include "ActsFatras/EventData/Barcode.hpp" @@ -33,12 +34,14 @@ #include using Acts::VectorHelpers::eta; +using Acts::VectorHelpers::perp; using Acts::VectorHelpers::phi; using Acts::VectorHelpers::theta; -ActsExamples::RootTrackParameterWriter::RootTrackParameterWriter( - const ActsExamples::RootTrackParameterWriter::Config& config, - Acts::Logging::Level level) +namespace ActsExamples { + +RootTrackParameterWriter::RootTrackParameterWriter( + const RootTrackParameterWriter::Config& config, Acts::Logging::Level level) : TrackParameterWriter(config.inputTrackParameters, "RootTrackParameterWriter", level), m_cfg(config) { @@ -83,37 +86,78 @@ ActsExamples::RootTrackParameterWriter::RootTrackParameterWriter( m_outputTree = new TTree(m_cfg.treeName.c_str(), m_cfg.treeName.c_str()); if (m_outputTree == nullptr) { throw std::bad_alloc(); - } else { - // The estimated track parameters - m_outputTree->Branch("event_nr", &m_eventNr); - m_outputTree->Branch("loc0", &m_loc0); - m_outputTree->Branch("loc1", &m_loc1); - m_outputTree->Branch("phi", &m_phi); - m_outputTree->Branch("theta", &m_theta); - m_outputTree->Branch("qop", &m_qop); - m_outputTree->Branch("time", &m_time); - m_outputTree->Branch("p", &m_p); - m_outputTree->Branch("pt", &m_pt); - m_outputTree->Branch("eta", &m_eta); - // The truth track parameters - m_outputTree->Branch("eventNr", &m_eventNr); - m_outputTree->Branch("t_loc0", &m_t_loc0); - m_outputTree->Branch("t_loc1", &m_t_loc1); - m_outputTree->Branch("t_phi", &m_t_phi); - m_outputTree->Branch("t_theta", &m_t_theta); - m_outputTree->Branch("t_qop", &m_t_qop); - m_outputTree->Branch("t_time", &m_t_time); - m_outputTree->Branch("truthMatched", &m_truthMatched); } + + m_outputTree->Branch("event_nr", &m_eventNr); + + m_outputTree->Branch("volumeId", &m_volumeId); + m_outputTree->Branch("layerId", &m_layerId); + m_outputTree->Branch("surfaceId", &m_surfaceId); + + // The estimated track parameters + m_outputTree->Branch("loc0", &m_loc0); + m_outputTree->Branch("loc1", &m_loc1); + m_outputTree->Branch("phi", &m_phi); + m_outputTree->Branch("theta", &m_theta); + m_outputTree->Branch("qop", &m_qop); + m_outputTree->Branch("time", &m_time); + + // The estimated track parameters errors + m_outputTree->Branch("err_loc0", &m_err_loc0); + m_outputTree->Branch("err_loc1", &m_err_loc1); + m_outputTree->Branch("err_phi", &m_err_phi); + m_outputTree->Branch("err_theta", &m_err_theta); + m_outputTree->Branch("err_qop", &m_err_qop); + m_outputTree->Branch("err_time", &m_err_time); + + // Some derived quantities from the estimated track parameters + m_outputTree->Branch("charge", &m_charge); + m_outputTree->Branch("p", &m_p); + m_outputTree->Branch("pt", &m_pt); + m_outputTree->Branch("eta", &m_eta); + + // The truth meta + m_outputTree->Branch("truthMatched", &m_t_matched); + m_outputTree->Branch("particleId", &m_t_particleId); + m_outputTree->Branch("nMajorityHits", &m_nMajorityHits); + + // The truth track parameters + m_outputTree->Branch("particleId", &m_t_particleId); + m_outputTree->Branch("t_loc0", &m_t_loc0); + m_outputTree->Branch("t_loc1", &m_t_loc1); + m_outputTree->Branch("t_phi", &m_t_phi); + m_outputTree->Branch("t_theta", &m_t_theta); + m_outputTree->Branch("t_qop", &m_t_qop); + m_outputTree->Branch("t_time", &m_t_time); + m_outputTree->Branch("t_charge", &m_t_charge); + m_outputTree->Branch("t_p", &m_t_p); + m_outputTree->Branch("t_pt", &m_t_pt); + m_outputTree->Branch("t_eta", &m_t_eta); + + // The residuals + m_outputTree->Branch("res_loc0", &m_res_loc0); + m_outputTree->Branch("res_loc1", &m_res_loc1); + m_outputTree->Branch("res_phi", &m_res_phi); + m_outputTree->Branch("res_theta", &m_res_theta); + m_outputTree->Branch("res_qop", &m_res_qop); + m_outputTree->Branch("res_time", &m_res_time); + + // The pulls + m_outputTree->Branch("pull_loc0", &m_pull_loc0); + m_outputTree->Branch("pull_loc1", &m_pull_loc1); + m_outputTree->Branch("pull_phi", &m_pull_phi); + m_outputTree->Branch("pull_theta", &m_pull_theta); + m_outputTree->Branch("pull_qop", &m_pull_qop); + m_outputTree->Branch("pull_time", &m_pull_time); } -ActsExamples::RootTrackParameterWriter::~RootTrackParameterWriter() { +RootTrackParameterWriter::~RootTrackParameterWriter() { if (m_outputFile != nullptr) { m_outputFile->Close(); } } -ActsExamples::ProcessCode ActsExamples::RootTrackParameterWriter::finalize() { +ProcessCode RootTrackParameterWriter::finalize() { m_outputFile->cd(); m_outputTree->Write(); m_outputFile->Close(); @@ -124,9 +168,8 @@ ActsExamples::ProcessCode ActsExamples::RootTrackParameterWriter::finalize() { return ProcessCode::SUCCESS; } -ActsExamples::ProcessCode ActsExamples::RootTrackParameterWriter::writeT( - const ActsExamples::AlgorithmContext& ctx, - const TrackParametersContainer& trackParams) { +ProcessCode RootTrackParameterWriter::writeT( + const AlgorithmContext& ctx, const TrackParametersContainer& trackParams) { // Read additional input collections const auto& protoTracks = m_inputProtoTracks(ctx); const auto& particles = m_inputParticles(ctx); @@ -134,6 +177,8 @@ ActsExamples::ProcessCode ActsExamples::RootTrackParameterWriter::writeT( const auto& hitParticlesMap = m_inputMeasurementParticlesMap(ctx); const auto& hitSimHitsMap = m_inputMeasurementSimHitsMap(ctx); + std::vector particleHitCounts; + // Exclusive access to the tree while writing std::lock_guard lock(m_writeMutex); @@ -144,65 +189,135 @@ ActsExamples::ProcessCode ActsExamples::RootTrackParameterWriter::writeT( // Loop over the estimated track parameters for (std::size_t iparams = 0; iparams < trackParams.size(); ++iparams) { + const auto& params = trackParams[iparams]; // The reference surface of the parameters, i.e. also the reference surface // of the first space point - const auto& surface = trackParams[iparams].referenceSurface(); - // The estimated bound parameters vector - const auto params = trackParams[iparams].parameters(); - m_loc0 = params[Acts::eBoundLoc0]; - m_loc1 = params[Acts::eBoundLoc1]; - m_phi = params[Acts::eBoundPhi]; - m_theta = params[Acts::eBoundTheta]; - m_qop = params[Acts::eBoundQOverP]; - m_time = params[Acts::eBoundTime]; - m_p = std::abs(1.0 / m_qop); - m_pt = m_p * std::sin(m_theta); - m_eta = std::atanh(std::cos(m_theta)); + const auto& surface = params.referenceSurface(); + + m_volumeId = surface.geometryId().volume(); + m_layerId = surface.geometryId().layer(); + m_surfaceId = surface.geometryId().sensitive(); + + m_loc0 = params.parameters()[Acts::eBoundLoc0]; + m_loc1 = params.parameters()[Acts::eBoundLoc1]; + m_phi = params.parameters()[Acts::eBoundPhi]; + m_theta = params.parameters()[Acts::eBoundTheta]; + m_qop = params.parameters()[Acts::eBoundQOverP]; + m_time = params.parameters()[Acts::eBoundTime]; + + auto getError = [¶ms](std::size_t idx) -> float { + if (!params.covariance().has_value()) { + return NaNfloat; + } + const auto& cov = *params.covariance(); + if (cov(idx, idx) < 0) { + return NaNfloat; + } + return std::sqrt(cov(idx, idx)); + }; + + m_err_loc0 = getError(Acts::eBoundLoc0); + m_err_loc1 = getError(Acts::eBoundLoc1); + m_err_phi = getError(Acts::eBoundPhi); + m_err_theta = getError(Acts::eBoundTheta); + m_err_qop = getError(Acts::eBoundQOverP); + m_err_time = getError(Acts::eBoundTime); + + m_charge = static_cast(params.charge()); + m_p = params.absoluteMomentum(); + m_pt = params.transverseMomentum(); + m_eta = eta(params.momentum()); // Get the proto track from which the track parameters are estimated const auto& ptrack = protoTracks[iparams]; - std::vector particleHitCounts; identifyContributingParticles(hitParticlesMap, ptrack, particleHitCounts); - m_truthMatched = false; + if (particleHitCounts.size() == 1) { - m_truthMatched = true; - } - // Get the index of the first space point - const auto& hitIdx = ptrack.front(); - // Get the sim hits via the measurement to sim hits map - auto indices = makeRange(hitSimHitsMap.equal_range(hitIdx)); - auto [truthLocal, truthPos4, truthUnitDir] = - averageSimHits(ctx.geoContext, surface, simHits, indices, logger()); - // Get the truth track parameter at the first space point - m_t_loc0 = truthLocal[Acts::ePos0]; - m_t_loc1 = truthLocal[Acts::ePos1]; - m_t_phi = phi(truthUnitDir); - m_t_theta = theta(truthUnitDir); - m_t_time = truthPos4[Acts::eTime]; - // momentum averaging makes even less sense than averaging position and - // direction. use the first momentum or set q/p to zero - if (!indices.empty()) { - // we assume that the indices are within valid ranges so we do not - // need to check their validity again. - const auto simHitIdx0 = indices.begin()->second; - const auto& simHit0 = *simHits.nth(simHitIdx0); - const auto p = - simHit0.momentum4Before().template segment<3>(Acts::eMom0).norm(); - const auto& particleId = simHit0.particleId(); - // The truth charge has to be retrieved from the sim particle - auto ip = particles.find(particleId); - if (ip != particles.end()) { - const auto& particle = *ip; - m_t_charge = static_cast(particle.charge()); - if (p != 0.0) { - m_t_qop = m_t_charge / p; + m_t_matched = true; + m_t_particleId = particleHitCounts.front().particleId.value(); + m_nMajorityHits = particleHitCounts.front().hitCount; + + // Get the index of the first space point + const auto& hitIdx = ptrack.front(); + // Get the sim hits via the measurement to sim hits map + auto indices = makeRange(hitSimHitsMap.equal_range(hitIdx)); + auto [truthLocal, truthPos4, truthUnitDir] = + averageSimHits(ctx.geoContext, surface, simHits, indices, logger()); + + // Get the truth track parameter at the first space point + m_t_loc0 = truthLocal[Acts::ePos0]; + m_t_loc1 = truthLocal[Acts::ePos1]; + m_t_phi = phi(truthUnitDir); + m_t_theta = theta(truthUnitDir); + m_t_qop = NaNfloat; + m_t_time = truthPos4[Acts::eTime]; + + m_t_charge = 0; + m_t_p = NaNfloat; + m_t_pt = NaNfloat; + m_t_eta = eta(truthUnitDir); + + // momentum averaging makes even less sense than averaging position and + // direction. use the first momentum or set q/p to zero + if (!indices.empty()) { + // we assume that the indices are within valid ranges so we do not + // need to check their validity again. + const auto simHitIdx0 = indices.begin()->second; + const auto& simHit0 = *simHits.nth(simHitIdx0); + const auto p = + simHit0.momentum4Before().template segment<3>(Acts::eMom0); + const auto& particleId = simHit0.particleId(); + // The truth charge has to be retrieved from the sim particle + if (auto ip = particles.find(particleId); ip != particles.end()) { + const auto& particle = *ip; + m_t_charge = static_cast(particle.charge()); + m_t_qop = particle.hypothesis().qOverP(p.norm(), particle.charge()); + m_t_p = p.norm(); + m_t_pt = perp(p); } else { - m_t_qop = std::numeric_limits::quiet_NaN(); + ACTS_WARNING("Truth particle with barcode " << particleId << "=" + << particleId.value() + << " not found!"); } - } else { - ACTS_DEBUG("Truth particle with barcode " - << particleId << "=" << particleId.value() << " not found!"); } + + m_res_loc0 = m_loc0 - m_t_loc0; + m_res_loc1 = m_loc1 - m_t_loc1; + m_res_phi = Acts::detail::difference_periodic( + m_phi, m_t_phi, static_cast(2 * M_PI)); + m_res_theta = m_theta - m_t_theta; + m_res_qop = m_qop - m_t_qop; + m_res_time = m_time - m_t_time; + + auto getPull = [](float res, float err) -> float { + if (std::isnan(err) || std::abs(err) < 1e-6) { + return NaNfloat; + } + return res / err; + }; + + m_pull_loc0 = getPull(m_res_loc0, m_err_loc0); + m_pull_loc1 = getPull(m_res_loc1, m_err_loc1); + m_pull_phi = getPull(m_res_phi, m_err_phi); + m_pull_theta = getPull(m_res_theta, m_err_theta); + m_pull_qop = getPull(m_res_qop, m_err_qop); + m_pull_time = getPull(m_res_time, m_err_time); + } else { + m_t_matched = false; + m_t_particleId = 0; + m_nMajorityHits = 0; + + m_t_loc0 = NaNfloat; + m_t_loc1 = NaNfloat; + m_t_phi = NaNfloat; + m_t_theta = NaNfloat; + m_t_qop = NaNfloat; + m_t_time = NaNfloat; + + m_t_charge = 0; + m_t_p = NaNfloat; + m_t_pt = NaNfloat; + m_t_eta = NaNfloat; } m_outputTree->Fill(); @@ -210,3 +325,5 @@ ActsExamples::ProcessCode ActsExamples::RootTrackParameterWriter::writeT( return ProcessCode::SUCCESS; } + +} // namespace ActsExamples diff --git a/Examples/Python/tests/root_file_hashes.txt b/Examples/Python/tests/root_file_hashes.txt index 7d86645f683..368a7c45ef3 100644 --- a/Examples/Python/tests/root_file_hashes.txt +++ b/Examples/Python/tests/root_file_hashes.txt @@ -3,16 +3,16 @@ test_fatras__particles_simulation.root: bc970873fef0c2efd86ed5413623802353d2cd04 test_fatras__hits.root: 6e4beb045fa1712c4d14c280ba33c3fa13e4aff9de88d55c3e32f62ad226f724 test_geant4__particles_simulation.root: 3c9c6265183b04c9d62ed5f2d0709c6dd74e33fbb53ac0aeb3274f6600257fc1 test_geant4__hits.root: adf5dcdf000a580412dc5089e17460897d6535c978eafa021584ba4281d0a1ac -test_seeding__estimatedparams.root: d18718704a9ba5bf6627a0e4bf60581c640d80d66fc46dbd3e7c1d051ad99f71 +test_seeding__estimatedparams.root: ff9b876f07de7e352c71f11dccdbaed93b6ef0328bba40c4f6c9fdb2090f386e test_seeding__performance_seeding.root: 992f9c611d30dde0d3f3ab676bab19ada61ab6a4442828e27b65ec5e5b7a2880 test_seeding__particles.root: 7855b021f39ad238bca098e4282667be0666f2d1630e5bcb9d51d3b5ee39fa14 test_seeding__particles_simulation.root: 87d9c6c82422ca381a17735096b36c547eacb5dda2f26d7377614bd97a70ab1a -test_hashing_seeding__estimatedparams.root: fa3a1af46651ea5982dfbf1a34dc81770e185a4b0be8e7a0666075987881d7fb -test_seeding_orthogonal__estimatedparams.root: 974b9f45c43b81504c24439b47895a5d0a0501f205d2c377923d303fc10adc31 +test_hashing_seeding__estimatedparams.root: 77658e952bbceb5a61b0235c76f4dbb7bb5994a8ef32b769197ddcfe46ade241 +test_seeding_orthogonal__estimatedparams.root: 47718b57d59d9b9686bd7a7bb5958d668c0a953ad6b094601ca7b27f530db786 test_seeding_orthogonal__performance_seeding.root: 60fbedcf5cb2b37cd8e526251940564432890d3a159d231ed819e915a904682c test_seeding_orthogonal__particles.root: 7855b021f39ad238bca098e4282667be0666f2d1630e5bcb9d51d3b5ee39fa14 test_seeding_orthogonal__particles_simulation.root: 87d9c6c82422ca381a17735096b36c547eacb5dda2f26d7377614bd97a70ab1a -test_itk_seeding__estimatedparams.root: 7e3acb54fcaabae19fe5b8601bd5bd3e5f65d6a1e6d5a2b59c4e8affd27e526c +test_itk_seeding__estimatedparams.root: 1266e99f1885a3c74a58471a4d6e750ad8fd87c7372a9ab929e8b99603b42ac8 test_itk_seeding__performance_seeding.root: 78ebda54cd0f026ba4b7f316724ffd946de56a932735914baf1b7bba9505c29d test_itk_seeding__particles.root: 0b6f4ad438010ac48803d48ed98e80b5e87d310bae6c2c02b16cd94d7a4d7d07 test_itk_seeding__particles_simulation.root: ef0246069aa697019f28a8b270a68de95312cae5f2f2c74848566c3ce4f70363