diff --git a/JugAlgo/src/components/AlgoServiceSvc.cpp b/JugAlgo/src/components/AlgoServiceSvc.cpp index 4af79335ad335f3c5018b87e2576becb6d6e3a15..90cb775c357133e40b8961632947dee5fb9b3d1d 100644 --- a/JugAlgo/src/components/AlgoServiceSvc.cpp +++ b/JugAlgo/src/components/AlgoServiceSvc.cpp @@ -9,6 +9,7 @@ #include #include +#include #include // too many services? :P @@ -22,6 +23,7 @@ public: private: SmartIF m_geoSvc; + Gaudi::Property m_randomSeed{this, "randomSeed", 1}; }; DECLARE_COMPONENT(AlgoServiceSvc) @@ -83,6 +85,14 @@ StatusCode AlgoServiceSvc::initialize() { info() << "Setting up algorithms::GeoSvc" << endmsg; auto* geo = static_cast(svc); geo->init(m_geoSvc->detector()); + } else if (name == algorithms::RandomSvc::kName) { + // setup random service + info() << "Setting up algorithms::RandomSvc\n" + << " --> using internal STL 64-bit MT engine\n" + << " --> seed set to" << m_randomSeed << endmsg; + auto* rnd = static_cast(svc); + rnd->setProperty("seed", m_randomSeed); + rnd->init(); } else { fatal() << "Unknown service encountered, please implement the necessary framework hooks" << endmsg; diff --git a/JugFast/CMakeLists.txt b/JugFast/CMakeLists.txt index 91b9c9b99413d59344ccde8d56d087329b4cf229..3b937ef3fef024a206e3f55e96415a405a32c373 100644 --- a/JugFast/CMakeLists.txt +++ b/JugFast/CMakeLists.txt @@ -12,6 +12,8 @@ gaudi_add_module(JugFastPlugins LINK Gaudi::GaudiAlgLib Gaudi::GaudiKernel JugBase + JugBase JugAlgo + algocore algotruth ROOT::Core ROOT::RIO ROOT::Tree EDM4HEP::edm4hep EDM4EIC::edm4eic diff --git a/JugFast/src/components/MC2SmearedParticle.cpp b/JugFast/src/components/MC2SmearedParticle.cpp index 6968e37fde8b36293b063a13ca3e5dd65e4f543f..148c870cc5cb0254b6f4ed07996a38ea5891c6fe 100644 --- a/JugFast/src/components/MC2SmearedParticle.cpp +++ b/JugFast/src/components/MC2SmearedParticle.cpp @@ -1,99 +1,8 @@ // SPDX-License-Identifier: LGPL-3.0-or-later // Copyright (C) 2022 Sylvester Joosten, Whitney Armstrong, Wouter Deconinck -#include "Gaudi/Algorithm.h" -#include "GaudiAlg/GaudiTool.h" -#include "GaudiAlg/Producer.h" -#include "GaudiAlg/Transformer.h" -#include "GaudiKernel/RndmGenerators.h" -#include -#include - -#include "JugBase/DataHandle.h" - -// Event Model related classes -#include "edm4hep/MCParticleCollection.h" -#include "edm4eic/ReconstructedParticleCollection.h" - -namespace Jug::Fast { - -class MC2SmearedParticle : public GaudiAlgorithm { -private: - DataHandle m_inputMCParticles{"MCParticles", Gaudi::DataHandle::Reader, this}; - DataHandle m_outputParticles{"SmearedReconstructedParticles", - Gaudi::DataHandle::Writer, this}; - Rndm::Numbers m_gaussDist; - Gaudi::Property m_smearing{this, "smearing", 0.01 /* 1 percent*/}; - -public: - MC2SmearedParticle(const std::string& name, ISvcLocator* svcLoc) : GaudiAlgorithm(name, svcLoc) { - declareProperty("inputParticles", m_inputMCParticles, "MCParticles"); - declareProperty("outputParticles", m_outputParticles, "SmearedReconstructedParticles"); - } - StatusCode initialize() override { - if (GaudiAlgorithm::initialize().isFailure()) { - return StatusCode::FAILURE; - } - - IRndmGenSvc* randSvc = svc("RndmGenSvc", true); - StatusCode sc = m_gaussDist.initialize(randSvc, Rndm::Gauss(1.0, m_smearing.value())); - if (!sc.isSuccess()) { - return StatusCode::FAILURE; - } - return StatusCode::SUCCESS; - } - StatusCode execute() override { - // input collection - const auto* const parts = m_inputMCParticles.get(); - // output collection - auto& out_parts = *(m_outputParticles.createAndPut()); - for (const auto& p : *parts) { - if (p.getGeneratorStatus() > 1) { - if (msgLevel(MSG::DEBUG)) { - debug() << "ignoring particle with generatorStatus = " << p.getGeneratorStatus() << endmsg; - } - continue; - } - - // for now just use total momentum smearing as this is the largest effect, - // ideally we should also smear the angles but this should be good enough - // for now. - const auto pvec = p.getMomentum(); - const auto pgen = std::hypot(pvec.x, pvec.y, pvec.z); - const auto momentum = pgen * m_gaussDist(); - // make sure we keep energy consistent - using MomType = decltype(edm4eic::ReconstructedParticle().getMomentum().x); - const MomType energy = std::sqrt(p.getEnergy() * p.getEnergy() - pgen * pgen + momentum * momentum); - const MomType px = p.getMomentum().x * momentum / pgen; - const MomType py = p.getMomentum().y * momentum / pgen; - const MomType pz = p.getMomentum().z * momentum / pgen; - - const MomType dpx = m_smearing.value() * px; - const MomType dpy = m_smearing.value() * py; - const MomType dpz = m_smearing.value() * pz; - const MomType dE = m_smearing.value() * energy; - // ignore covariance for now - // @TODO: vertex smearing - const MomType vx = p.getVertex().x; - const MomType vy = p.getVertex().y; - const MomType vz = p.getVertex().z; - - auto rec_part = out_parts.create(); - rec_part.setType(-1); // @TODO: determine type codes - rec_part.setEnergy(energy); - rec_part.setMomentum({px, py, pz}); - rec_part.setReferencePoint({vx, vy, vz}); // @FIXME: probably not what we want? - rec_part.setCharge(p.getCharge()); - rec_part.setMass(p.getMass()); - rec_part.setGoodnessOfPID(1); // Perfect PID - rec_part.setCovMatrix({dpx, dpy, dpz, dE}); - rec_part.setPDG(p.getPDG()); - } - return StatusCode::SUCCESS; - } -}; +#include +#include // NOLINTNEXTLINE(cppcoreguidelines-avoid-non-const-global-variables) -DECLARE_COMPONENT(MC2SmearedParticle) - -} // namespace Jug::Fast +JUGALGO_DEFINE_ALGORITHM(MC2SmearedParticle, algorithms::truth::MC2SmearedParticle, Jug::Fast) diff --git a/JugFast/src/components/ParticlesWithTruthPID.cpp b/JugFast/src/components/ParticlesWithTruthPID.cpp index 40177ae4900195466dc30afa243183977bcf9a01..face0d0728af278d6b21b936022caeea52b14f85 100644 --- a/JugFast/src/components/ParticlesWithTruthPID.cpp +++ b/JugFast/src/components/ParticlesWithTruthPID.cpp @@ -1,163 +1,8 @@ // SPDX-License-Identifier: LGPL-3.0-or-later // Copyright (C) 2022 Sylvester Joosten, Wouter Deconinck -#include -#include - -#include - -#include "Gaudi/Algorithm.h" -#include "GaudiAlg/GaudiTool.h" -#include "GaudiAlg/Producer.h" -#include "GaudiAlg/Transformer.h" -#include "GaudiKernel/RndmGenerators.h" - -#include "JugBase/DataHandle.h" - -// Event Model related classes -#include "edm4hep/MCParticleCollection.h" -#include "edm4eic/MCRecoParticleAssociationCollection.h" -#include "edm4eic/ReconstructedParticleCollection.h" -#include "edm4eic/TrackParametersCollection.h" -#include "edm4eic/vector_utils.h" - -namespace Jug::Fast { - -class ParticlesWithTruthPID : public GaudiAlgorithm { -private: - DataHandle m_inputTruthCollection{"inputMCParticles", Gaudi::DataHandle::Reader, this}; - DataHandle m_inputTrackCollection{"inputTrackParameters", Gaudi::DataHandle::Reader, - this}; - DataHandle m_outputParticleCollection{"ReconstructedParticles", - Gaudi::DataHandle::Writer, this}; - DataHandle m_outputAssocCollection{"MCRecoParticleAssociation", - Gaudi::DataHandle::Writer, this}; - - // Matching momentum tolerance requires 10% by default; - Gaudi::Property m_pRelativeTolerance{this, "pRelativeTolerance", {0.1}}; - // Matching phi tolerance of 10 mrad - Gaudi::Property m_phiTolerance{this, "phiTolerance", {0.030}}; - // Matchin eta tolerance of 0.1 - Gaudi::Property m_etaTolerance{this, "etaTolerance", {0.2}}; - -public: - ParticlesWithTruthPID(const std::string& name, ISvcLocator* svcLoc) : GaudiAlgorithm(name, svcLoc) { - declareProperty("inputMCParticles", m_inputTruthCollection, "MCParticles"); - declareProperty("inputTrackParameters", m_inputTrackCollection, "outputTrackParameters"); - declareProperty("outputParticles", m_outputParticleCollection, "ReconstructedParticles"); - declareProperty("outputAssociations", m_outputAssocCollection, "MCRecoParticleAssociation"); - } - StatusCode initialize() override { - if (GaudiAlgorithm::initialize().isFailure()) { - return StatusCode::FAILURE; - } - return StatusCode::SUCCESS; - } - StatusCode execute() override { - // input collection - const auto& mc = *(m_inputTruthCollection.get()); - const auto& tracks = *(m_inputTrackCollection.get()); - auto& part = *(m_outputParticleCollection.createAndPut()); - auto& assoc = *(m_outputAssocCollection.createAndPut()); - - const double sinPhiOver2Tolerance = sin(0.5 * m_phiTolerance); - std::vector consumed(mc.size(), false); - for (const auto& trk : tracks) { - const auto mom = edm4eic::sphericalToVector(1.0 / std::abs(trk.getQOverP()), trk.getTheta(), trk.getPhi()); - const auto charge_rec = trk.getCharge(); - // utility variables for matching - int best_match = -1; - double best_delta = std::numeric_limits::max(); - for (size_t ip = 0; ip < mc.size(); ++ip) { - const auto& mcpart = mc[ip]; - if (consumed[ip] || mcpart.getGeneratorStatus() > 1 || mcpart.getCharge() == 0 || - mcpart.getCharge() * charge_rec < 0) { - if (msgLevel(MSG::DEBUG)) { - debug() << "ignoring non-primary/neutral/opposite charge particle" << endmsg; - } - continue; - } - const auto& p = mcpart.getMomentum(); - const auto p_mag = std::hypot(p.x, p.y, p.z); - const auto p_phi = std::atan2(p.y, p.x); - const auto p_eta = std::atanh(p.z / p_mag); - const double dp_rel = std::abs((edm4eic::magnitude(mom) - p_mag) / p_mag); - // check the tolerance for sin(dphi/2) to avoid the hemisphere problem and allow - // for phi rollovers - const double dsphi = std::abs(sin(0.5 * (edm4eic::angleAzimuthal(mom) - p_phi))); - const double deta = std::abs((edm4eic::eta(mom) - p_eta)); - - if (dp_rel < m_pRelativeTolerance && deta < m_etaTolerance && dsphi < sinPhiOver2Tolerance) { - const double delta = - std::hypot(dp_rel / m_pRelativeTolerance, deta / m_etaTolerance, dsphi / sinPhiOver2Tolerance); - if (delta < best_delta) { - best_match = ip; - best_delta = delta; - } - } - } - auto rec_part = part.create(); - int32_t best_pid = 0; - auto referencePoint = rec_part.referencePoint(); - // float time = 0; - float mass = 0; - if (best_match >= 0) { - consumed[best_match] = true; - const auto& mcpart = mc[best_match]; - best_pid = mcpart.getPDG(); - referencePoint = { - static_cast(mcpart.getVertex().x), static_cast(mcpart.getVertex().y), - static_cast(mcpart.getVertex().z)}; // @TODO: not sure if vertex/reference poitn makes sense here - // time = mcpart.getTime(); - mass = mcpart.getMass(); - } - rec_part.setType(static_cast(best_match >= 0 ? 0 : -1)); // @TODO: determine type codes - rec_part.setEnergy(std::hypot(edm4eic::magnitude(mom), mass)); - rec_part.setMomentum(mom); - rec_part.setReferencePoint(referencePoint); - rec_part.setCharge(charge_rec); - rec_part.setMass(mass); - rec_part.setGoodnessOfPID(1); // perfect PID - rec_part.setPDG(best_pid); - // rec_part.covMatrix() // @TODO: covariance matrix on 4-momentum - // Also write MC <--> truth particle association if match was found - if (best_match >= 0) { - auto rec_assoc = assoc.create(); - rec_assoc.setRecID(rec_part.getObjectID().index); - rec_assoc.setSimID(mc[best_match].getObjectID().index); - rec_assoc.setWeight(1); - rec_assoc.setRec(rec_part); - //rec_assoc.setSim(mc[best_match]); - } - if (msgLevel(MSG::DEBUG)) { - if (best_match > 0) { - const auto& mcpart = mc[best_match]; - debug() << fmt::format("Matched track with MC particle {}\n", best_match) << endmsg; - debug() << fmt::format(" - Track: (mom: {}, theta: {}, phi: {}, charge: {})", edm4eic::magnitude(mom), - edm4eic::anglePolar(mom), edm4eic::angleAzimuthal(mom), charge_rec) - << endmsg; - const auto& p = mcpart.getMomentum(); - const auto p_mag = edm4eic::magnitude(p); - const auto p_phi = edm4eic::angleAzimuthal(p); - const auto p_theta = edm4eic::anglePolar(p); - debug() << fmt::format(" - MC particle: (mom: {}, theta: {}, phi: {}, charge: {}, type: {}", p_mag, p_theta, - p_phi, mcpart.getCharge(), mcpart.getPDG()) - << endmsg; - } else { - debug() << fmt::format("Did not find a good match for track \n") << endmsg; - debug() << fmt::format(" - Track: (mom: {}, theta: {}, phi: {}, charge: {})", edm4eic::magnitude(mom), - edm4eic::anglePolar(mom), edm4eic::angleAzimuthal(mom), charge_rec) - << endmsg; - } - } - } - - return StatusCode::SUCCESS; - } -}; // namespace Jug::Fast +#include +#include // NOLINTNEXTLINE(cppcoreguidelines-avoid-non-const-global-variables) -DECLARE_COMPONENT(ParticlesWithTruthPID) - -} // namespace Jug::Fast - +JUGALGO_DEFINE_ALGORITHM(ParticlesWithTruthPID, algorithms::truth::ParticlesWithTruthPID, Jug::Fast) diff --git a/external/algorithms/CMakeLists.txt b/external/algorithms/CMakeLists.txt index 393a3a0209b45e21e5da06c239d6233173ee04c4..68401fb774eea7e9e6272fdfebd2e88bc6051775 100644 --- a/external/algorithms/CMakeLists.txt +++ b/external/algorithms/CMakeLists.txt @@ -68,5 +68,5 @@ add_subdirectory(calorimetry) #add_subdirectory(far_forward) #add_subdirectory(pid) #add_subdirectory(tracking) -#add_subdirectory(truth) +add_subdirectory(truth) #add_subdirectory(utility) diff --git a/external/algorithms/calorimetry/include/algorithms/calorimetry/ClusterRecoCoG.h b/external/algorithms/calorimetry/include/algorithms/calorimetry/ClusterRecoCoG.h index b43090ac6a74ac32c783e372c97502669e488368..96502e07c0179676046ac30c0c8be6c155a9b454 100644 --- a/external/algorithms/calorimetry/include/algorithms/calorimetry/ClusterRecoCoG.h +++ b/external/algorithms/calorimetry/include/algorithms/calorimetry/ClusterRecoCoG.h @@ -10,7 +10,6 @@ #include #include -#include // Data types #include @@ -34,8 +33,6 @@ using ClusteringAlgorithm = Algorithm< */ class ClusterRecoCoG : public ClusteringAlgorithm { public: - using Input = ClusteringAlgorithm::Input; - using Output = ClusteringAlgorithm::Output; using WeightFunc = std::function; // TODO: get rid of "Collection" in names @@ -45,7 +42,7 @@ public: {"outputClusterCollection", "outputAssociations"}} {} void init(); - void process(const Input&, const Output&); + void process(const Input&, const Output&) const; private: edm4eic::MutableCluster reconstruct(const edm4eic::ProtoCluster&) const; diff --git a/external/algorithms/calorimetry/src/ClusterRecoCoG.cpp b/external/algorithms/calorimetry/src/ClusterRecoCoG.cpp index 95018716124fc4c2351602bfa754745eb6153e76..a456cd3f43ed9542cef710c852357ecbbedf9712 100644 --- a/external/algorithms/calorimetry/src/ClusterRecoCoG.cpp +++ b/external/algorithms/calorimetry/src/ClusterRecoCoG.cpp @@ -54,7 +54,7 @@ void ClusterRecoCoG::init() { } void ClusterRecoCoG::process(const ClusterRecoCoG::Input& input, - const ClusterRecoCoG::Output& output) { + const ClusterRecoCoG::Output& output) const { const auto [proto, opt_simhits] = input; auto [clusters, opt_assoc] = output; diff --git a/external/algorithms/core/include/algorithms/detail/random.h b/external/algorithms/core/include/algorithms/detail/random.h new file mode 100644 index 0000000000000000000000000000000000000000..cc43d7545b8e63e8734adb550f1bf4747d3c3237 --- /dev/null +++ b/external/algorithms/core/include/algorithms/detail/random.h @@ -0,0 +1,48 @@ +#pragma once + +#include +#include +#include + +namespace algorithms::detail { + +// Auto-refreshing cached sequence from an underlying random engine allowing for multiple instances +// to be evaluated in parallel. Specs: +// - GenFunc is required to return N random numbers between 0 and +// std::numeric_limits::max() +// - GenFunc is responsible to deal with possible simultaneous access by multiple +// instances of CachedBitGenerator (required to be thread-safe). +// - The owner of a CachedGenerator instance is responsible to prevent parallel +// calls to ::operator() +// Implements the uniform_random_bit_generator concept +class CachedBitGenerator { +public: + using value_type = uint_fast64_t; + using GenFunc = std::function(size_t /* N */)>; + CachedBitGenerator(const GenFunc& gen, const size_t cache_size) + // index starts at the end of the (empty) cache to force an immediate refresh + // on first access + : m_gen{gen}, m_cache(cache_size), m_index{cache_size} {} + + value_type operator()() { + if (m_index >= m_cache.size()) { + refresh(); + } + return m_cache[m_index++]; + } + + constexpr value_type min() const { return 0; } + constexpr value_type max() const { return std::numeric_limits::max(); } + +private: + void refresh() { + m_cache = m_gen(m_cache.size()); + m_index = 0; + } + + GenFunc m_gen; + std::vector m_cache; + size_t m_index; +}; + +} // namespace algorithms::detail diff --git a/external/algorithms/core/include/algorithms/random.h b/external/algorithms/core/include/algorithms/random.h new file mode 100644 index 0000000000000000000000000000000000000000..6826aada45e2c9759b3688417c616f18191f8c68 --- /dev/null +++ b/external/algorithms/core/include/algorithms/random.h @@ -0,0 +1,115 @@ +#pragma once + +#include +#include +#include +#include + +#include +#include +#include + +namespace algorithms { + +// Random Engine callback function: +// - Signature: std::function(size_t N)> --> generates a vector +// of N numbers +// - RandomEngineCB is required to return N random numbers between 0 and +// std::numeric_limits::max() +// - RandomEngineCB is responsible to deal with possible simultaneous access by multiple +// Generator instances (required to be thread-safe). +using RandomEngineCB = detail::CachedBitGenerator::GenFunc; + +// thread-safe generator front-end. Requires that the underlying random engine used by +// the RandomSvc is thread-safe. +class Generator { +public: + Generator(const RandomEngineCB& gen, const size_t cache_size) : m_gen{gen, cache_size} {} + + template Int uniform_int(const Int min, const Int max) const { + std::uniform_int_distribution d{min, max}; + std::lock_guard lock{m_mutex}; + return d(m_gen); + } + template Float uniform_double(const Float min, const Float max) const { + std::uniform_real_distribution d{min, max}; + std::lock_guard lock{m_mutex}; + return d(m_gen); + } + template Int poisson(const Int mean) const { + std::poisson_distribution d{mean}; + std::lock_guard lock{m_mutex}; + return d(m_gen); + } + template Float exponential(const Float lambda) const { + std::exponential_distribution d{lambda}; + std::lock_guard lock{m_mutex}; + return d(m_gen); + } + template Float gaussian(const Float mu, const Float sigma) const { + std::normal_distribution d{mu, sigma}; + std::lock_guard lock{m_mutex}; + return d(m_gen); + } + +private: + mutable detail::CachedBitGenerator m_gen; + mutable std::mutex m_mutex; +}; + +// Random service that creates multiple Generators that are linked to a single random +// engine. The Generators are safe to be used in parallel as long as the Engine itself is +// thread-safe (this is a hard requirement for MT). The Generators avoid unnecesary locking +// by running off a auto-refreshing cached random sequence. +class RandomSvc : public LoggedService { +public: + using value_type = detail::CachedBitGenerator::value_type; + + Generator generator() { return {m_gen, m_cache_size}; } +// FIXME fix the CMake setup so these are properly found in Gaudi +#if 0 + void init(); + void init(const RandomEngineCB& gen); +#endif + void init() { + if (m_seed.hasValue()) { + info() << "Custom random seed requested: " << m_seed << endmsg; + m_gen = createEngine(m_seed); + } + } + void init(const RandomEngineCB& gen) { + init(); + info() << "Loading external generator function." << endmsg; + m_gen = gen; + if (m_seed.hasValue()) { + warning() << "Custom random seed request ignored when using external generator function" + << endmsg; + } + } + +#if 0 +private: + RandomEngineCB createEngine(const size_t seed = 1); +#endif + RandomEngineCB createEngine(const size_t seed = 1) { + return [=](const size_t size) { + static std::mutex m; + static std::mt19937_64 gen{seed}; + std::lock_guard lock{m}; + std::vector ret(size); + std::generate(ret.begin(), ret.end(), gen); + return ret; + }; + } + // end of FIXME + +private: + RandomEngineCB m_gen{createEngine()}; + Property m_seed{this, "seed"}; + Property m_cache_size{this, "cacheSize", 1024}; + std::mutex m_mutex; + + ALGORITHMS_DEFINE_LOGGED_SERVICE(RandomSvc) +}; + +} // namespace algorithms diff --git a/external/algorithms/core/src/dummy.cpp b/external/algorithms/core/src/dummy.cpp index a264023dc3529dc6f4ed6effad7e953d9a7f33a5..5041865fb3ac09ad5741daee47bc65ede1043cf5 100644 --- a/external/algorithms/core/src/dummy.cpp +++ b/external/algorithms/core/src/dummy.cpp @@ -6,6 +6,7 @@ #include #include +#include #include diff --git a/external/algorithms/core/src/random.cpp b/external/algorithms/core/src/random.cpp new file mode 100644 index 0000000000000000000000000000000000000000..711e35abe35a4aa1ad2158c449db9bf68de0dda1 --- /dev/null +++ b/external/algorithms/core/src/random.cpp @@ -0,0 +1,33 @@ +#include + +namespace algorithms { + +// FIXME fix the CMake setup so these are properly found in Gaudi... +#if 0 +void RandomSvc::init() { + if (m_seed.hasValue()) { + info() << "Custom random seed requested: " << m_seed << endmsg; + m_gen = createEngine(m_seed); + } +} +void RandomSvc::init(const RandomEngineCB& gen) { + init(); + info() << "Loading external generator function." << endmsg; + m_gen = gen; + if (m_seed.hasValue()) { + warning() << "Custom random seed request ignored when using external generator function" + << endmsg; + } +} +RandomEngineCB RandomSvc::createEngine(const size_t seed) { + return [=](const size_t size) { + static std::mutex m; + static std::mt19937_64 gen{seed}; + std::lock_guard lock{m}; + std::vector ret{size}; + std::generate(ret.begin(), ret.end(), gen); + return ret; + }; +} +#endif +} // namespace algorithms diff --git a/external/algorithms/truth/CMakeLists.txt b/external/algorithms/truth/CMakeLists.txt index cca38a857bdf4643807974bdc2b334e728c88f90..7c465f9298fa665bd4409bcf1062e0c2648f2f3b 100644 --- a/external/algorithms/truth/CMakeLists.txt +++ b/external/algorithms/truth/CMakeLists.txt @@ -12,6 +12,8 @@ set(TARGETS ${TARGETS} ${LIBRARY} PARENT_SCOPE) # FIXME: adding one by one #file(GLOB SRC CONFIGURE_DEPENDS src/*.cpp) set(SRC + src/MC2SmearedParticle.cpp + src/ParticlesWithTruthPID.cpp ) add_library(${LIBRARY} SHARED ${SRC}) diff --git a/external/algorithms/truth/include/algorithms/truth/MC2SmearedParticle.cpp b/external/algorithms/truth/include/algorithms/truth/MC2SmearedParticle.cpp deleted file mode 100644 index 6968e37fde8b36293b063a13ca3e5dd65e4f543f..0000000000000000000000000000000000000000 --- a/external/algorithms/truth/include/algorithms/truth/MC2SmearedParticle.cpp +++ /dev/null @@ -1,99 +0,0 @@ -// SPDX-License-Identifier: LGPL-3.0-or-later -// Copyright (C) 2022 Sylvester Joosten, Whitney Armstrong, Wouter Deconinck - -#include "Gaudi/Algorithm.h" -#include "GaudiAlg/GaudiTool.h" -#include "GaudiAlg/Producer.h" -#include "GaudiAlg/Transformer.h" -#include "GaudiKernel/RndmGenerators.h" -#include -#include - -#include "JugBase/DataHandle.h" - -// Event Model related classes -#include "edm4hep/MCParticleCollection.h" -#include "edm4eic/ReconstructedParticleCollection.h" - -namespace Jug::Fast { - -class MC2SmearedParticle : public GaudiAlgorithm { -private: - DataHandle m_inputMCParticles{"MCParticles", Gaudi::DataHandle::Reader, this}; - DataHandle m_outputParticles{"SmearedReconstructedParticles", - Gaudi::DataHandle::Writer, this}; - Rndm::Numbers m_gaussDist; - Gaudi::Property m_smearing{this, "smearing", 0.01 /* 1 percent*/}; - -public: - MC2SmearedParticle(const std::string& name, ISvcLocator* svcLoc) : GaudiAlgorithm(name, svcLoc) { - declareProperty("inputParticles", m_inputMCParticles, "MCParticles"); - declareProperty("outputParticles", m_outputParticles, "SmearedReconstructedParticles"); - } - StatusCode initialize() override { - if (GaudiAlgorithm::initialize().isFailure()) { - return StatusCode::FAILURE; - } - - IRndmGenSvc* randSvc = svc("RndmGenSvc", true); - StatusCode sc = m_gaussDist.initialize(randSvc, Rndm::Gauss(1.0, m_smearing.value())); - if (!sc.isSuccess()) { - return StatusCode::FAILURE; - } - return StatusCode::SUCCESS; - } - StatusCode execute() override { - // input collection - const auto* const parts = m_inputMCParticles.get(); - // output collection - auto& out_parts = *(m_outputParticles.createAndPut()); - for (const auto& p : *parts) { - if (p.getGeneratorStatus() > 1) { - if (msgLevel(MSG::DEBUG)) { - debug() << "ignoring particle with generatorStatus = " << p.getGeneratorStatus() << endmsg; - } - continue; - } - - // for now just use total momentum smearing as this is the largest effect, - // ideally we should also smear the angles but this should be good enough - // for now. - const auto pvec = p.getMomentum(); - const auto pgen = std::hypot(pvec.x, pvec.y, pvec.z); - const auto momentum = pgen * m_gaussDist(); - // make sure we keep energy consistent - using MomType = decltype(edm4eic::ReconstructedParticle().getMomentum().x); - const MomType energy = std::sqrt(p.getEnergy() * p.getEnergy() - pgen * pgen + momentum * momentum); - const MomType px = p.getMomentum().x * momentum / pgen; - const MomType py = p.getMomentum().y * momentum / pgen; - const MomType pz = p.getMomentum().z * momentum / pgen; - - const MomType dpx = m_smearing.value() * px; - const MomType dpy = m_smearing.value() * py; - const MomType dpz = m_smearing.value() * pz; - const MomType dE = m_smearing.value() * energy; - // ignore covariance for now - // @TODO: vertex smearing - const MomType vx = p.getVertex().x; - const MomType vy = p.getVertex().y; - const MomType vz = p.getVertex().z; - - auto rec_part = out_parts.create(); - rec_part.setType(-1); // @TODO: determine type codes - rec_part.setEnergy(energy); - rec_part.setMomentum({px, py, pz}); - rec_part.setReferencePoint({vx, vy, vz}); // @FIXME: probably not what we want? - rec_part.setCharge(p.getCharge()); - rec_part.setMass(p.getMass()); - rec_part.setGoodnessOfPID(1); // Perfect PID - rec_part.setCovMatrix({dpx, dpy, dpz, dE}); - rec_part.setPDG(p.getPDG()); - } - return StatusCode::SUCCESS; - } -}; - -// NOLINTNEXTLINE(cppcoreguidelines-avoid-non-const-global-variables) -DECLARE_COMPONENT(MC2SmearedParticle) - -} // namespace Jug::Fast diff --git a/external/algorithms/truth/include/algorithms/truth/MC2SmearedParticle.h b/external/algorithms/truth/include/algorithms/truth/MC2SmearedParticle.h new file mode 100644 index 0000000000000000000000000000000000000000..e7ee0dfb2fdf64aaad0397a64584d49b3894e960 --- /dev/null +++ b/external/algorithms/truth/include/algorithms/truth/MC2SmearedParticle.h @@ -0,0 +1,31 @@ +// SPDX-License-Identifier: LGPL-3.0-or-later +// Copyright (C) 2022 Sylvester Joosten, Whitney Armstrong, Wouter Deconinck + +#include +#include + +// Event Model related classes +#include +#include + +namespace algorithms::truth { + +using MC2SmearedParticleAlgorithm = Algorithm, + Output>; + +class MC2SmearedParticle : public MC2SmearedParticleAlgorithm { +public: + MC2SmearedParticle(std::string_view name) + : MC2SmearedParticleAlgorithm{name, {"inputParticles"}, {"outputParticles"}} {} + + void init(); + void process(const Input&, const Output&) const; + +private: + Generator m_rng = RandomSvc::instance().generator(); + + Property m_smearing{this, "smearing", 0.01 /* 1 percent*/}; +}; + +} // namespace algorithms::truth + diff --git a/external/algorithms/truth/include/algorithms/truth/ParticlesWithTruthPID.cpp b/external/algorithms/truth/include/algorithms/truth/ParticlesWithTruthPID.cpp deleted file mode 100644 index 40177ae4900195466dc30afa243183977bcf9a01..0000000000000000000000000000000000000000 --- a/external/algorithms/truth/include/algorithms/truth/ParticlesWithTruthPID.cpp +++ /dev/null @@ -1,163 +0,0 @@ -// SPDX-License-Identifier: LGPL-3.0-or-later -// Copyright (C) 2022 Sylvester Joosten, Wouter Deconinck - -#include -#include - -#include - -#include "Gaudi/Algorithm.h" -#include "GaudiAlg/GaudiTool.h" -#include "GaudiAlg/Producer.h" -#include "GaudiAlg/Transformer.h" -#include "GaudiKernel/RndmGenerators.h" - -#include "JugBase/DataHandle.h" - -// Event Model related classes -#include "edm4hep/MCParticleCollection.h" -#include "edm4eic/MCRecoParticleAssociationCollection.h" -#include "edm4eic/ReconstructedParticleCollection.h" -#include "edm4eic/TrackParametersCollection.h" -#include "edm4eic/vector_utils.h" - -namespace Jug::Fast { - -class ParticlesWithTruthPID : public GaudiAlgorithm { -private: - DataHandle m_inputTruthCollection{"inputMCParticles", Gaudi::DataHandle::Reader, this}; - DataHandle m_inputTrackCollection{"inputTrackParameters", Gaudi::DataHandle::Reader, - this}; - DataHandle m_outputParticleCollection{"ReconstructedParticles", - Gaudi::DataHandle::Writer, this}; - DataHandle m_outputAssocCollection{"MCRecoParticleAssociation", - Gaudi::DataHandle::Writer, this}; - - // Matching momentum tolerance requires 10% by default; - Gaudi::Property m_pRelativeTolerance{this, "pRelativeTolerance", {0.1}}; - // Matching phi tolerance of 10 mrad - Gaudi::Property m_phiTolerance{this, "phiTolerance", {0.030}}; - // Matchin eta tolerance of 0.1 - Gaudi::Property m_etaTolerance{this, "etaTolerance", {0.2}}; - -public: - ParticlesWithTruthPID(const std::string& name, ISvcLocator* svcLoc) : GaudiAlgorithm(name, svcLoc) { - declareProperty("inputMCParticles", m_inputTruthCollection, "MCParticles"); - declareProperty("inputTrackParameters", m_inputTrackCollection, "outputTrackParameters"); - declareProperty("outputParticles", m_outputParticleCollection, "ReconstructedParticles"); - declareProperty("outputAssociations", m_outputAssocCollection, "MCRecoParticleAssociation"); - } - StatusCode initialize() override { - if (GaudiAlgorithm::initialize().isFailure()) { - return StatusCode::FAILURE; - } - return StatusCode::SUCCESS; - } - StatusCode execute() override { - // input collection - const auto& mc = *(m_inputTruthCollection.get()); - const auto& tracks = *(m_inputTrackCollection.get()); - auto& part = *(m_outputParticleCollection.createAndPut()); - auto& assoc = *(m_outputAssocCollection.createAndPut()); - - const double sinPhiOver2Tolerance = sin(0.5 * m_phiTolerance); - std::vector consumed(mc.size(), false); - for (const auto& trk : tracks) { - const auto mom = edm4eic::sphericalToVector(1.0 / std::abs(trk.getQOverP()), trk.getTheta(), trk.getPhi()); - const auto charge_rec = trk.getCharge(); - // utility variables for matching - int best_match = -1; - double best_delta = std::numeric_limits::max(); - for (size_t ip = 0; ip < mc.size(); ++ip) { - const auto& mcpart = mc[ip]; - if (consumed[ip] || mcpart.getGeneratorStatus() > 1 || mcpart.getCharge() == 0 || - mcpart.getCharge() * charge_rec < 0) { - if (msgLevel(MSG::DEBUG)) { - debug() << "ignoring non-primary/neutral/opposite charge particle" << endmsg; - } - continue; - } - const auto& p = mcpart.getMomentum(); - const auto p_mag = std::hypot(p.x, p.y, p.z); - const auto p_phi = std::atan2(p.y, p.x); - const auto p_eta = std::atanh(p.z / p_mag); - const double dp_rel = std::abs((edm4eic::magnitude(mom) - p_mag) / p_mag); - // check the tolerance for sin(dphi/2) to avoid the hemisphere problem and allow - // for phi rollovers - const double dsphi = std::abs(sin(0.5 * (edm4eic::angleAzimuthal(mom) - p_phi))); - const double deta = std::abs((edm4eic::eta(mom) - p_eta)); - - if (dp_rel < m_pRelativeTolerance && deta < m_etaTolerance && dsphi < sinPhiOver2Tolerance) { - const double delta = - std::hypot(dp_rel / m_pRelativeTolerance, deta / m_etaTolerance, dsphi / sinPhiOver2Tolerance); - if (delta < best_delta) { - best_match = ip; - best_delta = delta; - } - } - } - auto rec_part = part.create(); - int32_t best_pid = 0; - auto referencePoint = rec_part.referencePoint(); - // float time = 0; - float mass = 0; - if (best_match >= 0) { - consumed[best_match] = true; - const auto& mcpart = mc[best_match]; - best_pid = mcpart.getPDG(); - referencePoint = { - static_cast(mcpart.getVertex().x), static_cast(mcpart.getVertex().y), - static_cast(mcpart.getVertex().z)}; // @TODO: not sure if vertex/reference poitn makes sense here - // time = mcpart.getTime(); - mass = mcpart.getMass(); - } - rec_part.setType(static_cast(best_match >= 0 ? 0 : -1)); // @TODO: determine type codes - rec_part.setEnergy(std::hypot(edm4eic::magnitude(mom), mass)); - rec_part.setMomentum(mom); - rec_part.setReferencePoint(referencePoint); - rec_part.setCharge(charge_rec); - rec_part.setMass(mass); - rec_part.setGoodnessOfPID(1); // perfect PID - rec_part.setPDG(best_pid); - // rec_part.covMatrix() // @TODO: covariance matrix on 4-momentum - // Also write MC <--> truth particle association if match was found - if (best_match >= 0) { - auto rec_assoc = assoc.create(); - rec_assoc.setRecID(rec_part.getObjectID().index); - rec_assoc.setSimID(mc[best_match].getObjectID().index); - rec_assoc.setWeight(1); - rec_assoc.setRec(rec_part); - //rec_assoc.setSim(mc[best_match]); - } - if (msgLevel(MSG::DEBUG)) { - if (best_match > 0) { - const auto& mcpart = mc[best_match]; - debug() << fmt::format("Matched track with MC particle {}\n", best_match) << endmsg; - debug() << fmt::format(" - Track: (mom: {}, theta: {}, phi: {}, charge: {})", edm4eic::magnitude(mom), - edm4eic::anglePolar(mom), edm4eic::angleAzimuthal(mom), charge_rec) - << endmsg; - const auto& p = mcpart.getMomentum(); - const auto p_mag = edm4eic::magnitude(p); - const auto p_phi = edm4eic::angleAzimuthal(p); - const auto p_theta = edm4eic::anglePolar(p); - debug() << fmt::format(" - MC particle: (mom: {}, theta: {}, phi: {}, charge: {}, type: {}", p_mag, p_theta, - p_phi, mcpart.getCharge(), mcpart.getPDG()) - << endmsg; - } else { - debug() << fmt::format("Did not find a good match for track \n") << endmsg; - debug() << fmt::format(" - Track: (mom: {}, theta: {}, phi: {}, charge: {})", edm4eic::magnitude(mom), - edm4eic::anglePolar(mom), edm4eic::angleAzimuthal(mom), charge_rec) - << endmsg; - } - } - } - - return StatusCode::SUCCESS; - } -}; // namespace Jug::Fast - -// NOLINTNEXTLINE(cppcoreguidelines-avoid-non-const-global-variables) -DECLARE_COMPONENT(ParticlesWithTruthPID) - -} // namespace Jug::Fast - diff --git a/external/algorithms/truth/include/algorithms/truth/ParticlesWithTruthPID.h b/external/algorithms/truth/include/algorithms/truth/ParticlesWithTruthPID.h new file mode 100644 index 0000000000000000000000000000000000000000..d00fb76c202ad50fa58a5a5969a9713fbc4fc4b8 --- /dev/null +++ b/external/algorithms/truth/include/algorithms/truth/ParticlesWithTruthPID.h @@ -0,0 +1,38 @@ +// SPDX-License-Identifier: LGPL-3.0-or-later +// Copyright (C) 2022 Sylvester Joosten, Wouter Deconinck + +#include + +// Event Model related classes +#include +#include +#include +#include + +namespace algorithms::truth { + +using ParticlesWithTruthPIDAlgorithm = Algorithm< + Input, + Output>; + +class ParticlesWithTruthPID : public ParticlesWithTruthPIDAlgorithm { +public: + ParticlesWithTruthPID(std::string_view name) + : ParticlesWithTruthPIDAlgorithm{name, + {"inputMCParticles", "inputTrackParameters"}, + {"outputParticles", "outputAssociations"}} {} + + void init(); + void process(const Input&, const Output&) const; + +private: + // Matching momentum tolerance requires 10% by default; + Property m_pRelativeTolerance{this, "pRelativeTolerance", {0.1}}; + // Matching phi tolerance of 10 mrad + Property m_phiTolerance{this, "phiTolerance", {0.030}}; + // Matchin eta tolerance of 0.1 + Property m_etaTolerance{this, "etaTolerance", {0.2}}; +}; + +} // namespace algorithms::truth + diff --git a/external/algorithms/truth/src/MC2SmearedParticle.cpp b/external/algorithms/truth/src/MC2SmearedParticle.cpp index 6968e37fde8b36293b063a13ca3e5dd65e4f543f..89d9ccd3083d0a34bd6344a6764440aa3dbc56b7 100644 --- a/external/algorithms/truth/src/MC2SmearedParticle.cpp +++ b/external/algorithms/truth/src/MC2SmearedParticle.cpp @@ -1,99 +1,66 @@ // SPDX-License-Identifier: LGPL-3.0-or-later // Copyright (C) 2022 Sylvester Joosten, Whitney Armstrong, Wouter Deconinck -#include "Gaudi/Algorithm.h" -#include "GaudiAlg/GaudiTool.h" -#include "GaudiAlg/Producer.h" -#include "GaudiAlg/Transformer.h" -#include "GaudiKernel/RndmGenerators.h" -#include -#include - -#include "JugBase/DataHandle.h" +#include -// Event Model related classes -#include "edm4hep/MCParticleCollection.h" -#include "edm4eic/ReconstructedParticleCollection.h" +#include +#include -namespace Jug::Fast { +namespace algorithms::truth { -class MC2SmearedParticle : public GaudiAlgorithm { -private: - DataHandle m_inputMCParticles{"MCParticles", Gaudi::DataHandle::Reader, this}; - DataHandle m_outputParticles{"SmearedReconstructedParticles", - Gaudi::DataHandle::Writer, this}; - Rndm::Numbers m_gaussDist; - Gaudi::Property m_smearing{this, "smearing", 0.01 /* 1 percent*/}; +void MC2SmearedParticle::init() { + ; // do nothing +} -public: - MC2SmearedParticle(const std::string& name, ISvcLocator* svcLoc) : GaudiAlgorithm(name, svcLoc) { - declareProperty("inputParticles", m_inputMCParticles, "MCParticles"); - declareProperty("outputParticles", m_outputParticles, "SmearedReconstructedParticles"); - } - StatusCode initialize() override { - if (GaudiAlgorithm::initialize().isFailure()) { - return StatusCode::FAILURE; - } +void MC2SmearedParticle::process(const MC2SmearedParticle::Input& input, + const MC2SmearedParticle::Output& output) const { + const auto [parts] = input; + auto [out_parts] = output; - IRndmGenSvc* randSvc = svc("RndmGenSvc", true); - StatusCode sc = m_gaussDist.initialize(randSvc, Rndm::Gauss(1.0, m_smearing.value())); - if (!sc.isSuccess()) { - return StatusCode::FAILURE; - } - return StatusCode::SUCCESS; - } - StatusCode execute() override { - // input collection - const auto* const parts = m_inputMCParticles.get(); - // output collection - auto& out_parts = *(m_outputParticles.createAndPut()); - for (const auto& p : *parts) { - if (p.getGeneratorStatus() > 1) { - if (msgLevel(MSG::DEBUG)) { - debug() << "ignoring particle with generatorStatus = " << p.getGeneratorStatus() << endmsg; - } - continue; + for (const auto& p : *parts) { + if (p.getGeneratorStatus() > 1) { + if (aboveDebugThreshold()) { + debug() << "ignoring particle with generatorStatus = " << p.getGeneratorStatus() << endmsg; } + continue; + } - // for now just use total momentum smearing as this is the largest effect, - // ideally we should also smear the angles but this should be good enough - // for now. - const auto pvec = p.getMomentum(); - const auto pgen = std::hypot(pvec.x, pvec.y, pvec.z); - const auto momentum = pgen * m_gaussDist(); - // make sure we keep energy consistent - using MomType = decltype(edm4eic::ReconstructedParticle().getMomentum().x); - const MomType energy = std::sqrt(p.getEnergy() * p.getEnergy() - pgen * pgen + momentum * momentum); - const MomType px = p.getMomentum().x * momentum / pgen; - const MomType py = p.getMomentum().y * momentum / pgen; - const MomType pz = p.getMomentum().z * momentum / pgen; + // for now just use total momentum smearing as this is the largest effect, + // ideally we should also smear the angles but this should be good enough + // for now. + const auto pvec = p.getMomentum(); + const auto pgen = std::hypot(pvec.x, pvec.y, pvec.z); + const auto momentum = pgen * m_rng.gaussian(0., m_smearing); + // make sure we keep energy consistent + using MomType = decltype(edm4eic::ReconstructedParticle().getMomentum().x); + const MomType energy = + std::sqrt(p.getEnergy() * p.getEnergy() - pgen * pgen + momentum * momentum); + const MomType px = p.getMomentum().x * momentum / pgen; + const MomType py = p.getMomentum().y * momentum / pgen; + const MomType pz = p.getMomentum().z * momentum / pgen; - const MomType dpx = m_smearing.value() * px; - const MomType dpy = m_smearing.value() * py; - const MomType dpz = m_smearing.value() * pz; - const MomType dE = m_smearing.value() * energy; - // ignore covariance for now - // @TODO: vertex smearing - const MomType vx = p.getVertex().x; - const MomType vy = p.getVertex().y; - const MomType vz = p.getVertex().z; + const MomType dpx = m_smearing * px; + const MomType dpy = m_smearing * py; + const MomType dpz = m_smearing * pz; + const MomType dE = m_smearing * energy; + // ignore covariance for now + // @TODO: vertex smearing + const MomType vx = p.getVertex().x; + const MomType vy = p.getVertex().y; + const MomType vz = p.getVertex().z; - auto rec_part = out_parts.create(); - rec_part.setType(-1); // @TODO: determine type codes - rec_part.setEnergy(energy); - rec_part.setMomentum({px, py, pz}); - rec_part.setReferencePoint({vx, vy, vz}); // @FIXME: probably not what we want? - rec_part.setCharge(p.getCharge()); - rec_part.setMass(p.getMass()); - rec_part.setGoodnessOfPID(1); // Perfect PID - rec_part.setCovMatrix({dpx, dpy, dpz, dE}); - rec_part.setPDG(p.getPDG()); - } - return StatusCode::SUCCESS; + auto rec_part = out_parts->create(); + rec_part.setType(-1); // @TODO: determine type codes + rec_part.setEnergy(energy); + rec_part.setMomentum({px, py, pz}); + rec_part.setReferencePoint({vx, vy, vz}); // @FIXME: probably not what we want? + rec_part.setCharge(p.getCharge()); + rec_part.setMass(p.getMass()); + rec_part.setGoodnessOfPID(1); // Perfect PID + rec_part.setCovMatrix({dpx, dpy, dpz, dE}); + rec_part.setPDG(p.getPDG()); } -}; +} -// NOLINTNEXTLINE(cppcoreguidelines-avoid-non-const-global-variables) -DECLARE_COMPONENT(MC2SmearedParticle) +} // namespace algorithms::truth -} // namespace Jug::Fast diff --git a/external/algorithms/truth/src/ParticlesWithTruthPID.cpp b/external/algorithms/truth/src/ParticlesWithTruthPID.cpp index 40177ae4900195466dc30afa243183977bcf9a01..4cbe18e218faa23a2d112454bc3b3585160a497c 100644 --- a/external/algorithms/truth/src/ParticlesWithTruthPID.cpp +++ b/external/algorithms/truth/src/ParticlesWithTruthPID.cpp @@ -1,163 +1,126 @@ // SPDX-License-Identifier: LGPL-3.0-or-later // Copyright (C) 2022 Sylvester Joosten, Wouter Deconinck +#include + #include #include - #include -#include "Gaudi/Algorithm.h" -#include "GaudiAlg/GaudiTool.h" -#include "GaudiAlg/Producer.h" -#include "GaudiAlg/Transformer.h" -#include "GaudiKernel/RndmGenerators.h" - -#include "JugBase/DataHandle.h" +#include -// Event Model related classes -#include "edm4hep/MCParticleCollection.h" -#include "edm4eic/MCRecoParticleAssociationCollection.h" -#include "edm4eic/ReconstructedParticleCollection.h" -#include "edm4eic/TrackParametersCollection.h" -#include "edm4eic/vector_utils.h" +namespace algorithms::truth { -namespace Jug::Fast { +void ParticlesWithTruthPID::init() { + ; // do nothing +} +void ParticlesWithTruthPID::process(const ParticlesWithTruthPID::Input& input, + const ParticlesWithTruthPID::Output& output) const { + const auto [mc_ptr, tracks_ptr] = input; + auto [part_ptr, assoc_ptr] = output; -class ParticlesWithTruthPID : public GaudiAlgorithm { -private: - DataHandle m_inputTruthCollection{"inputMCParticles", Gaudi::DataHandle::Reader, this}; - DataHandle m_inputTrackCollection{"inputTrackParameters", Gaudi::DataHandle::Reader, - this}; - DataHandle m_outputParticleCollection{"ReconstructedParticles", - Gaudi::DataHandle::Writer, this}; - DataHandle m_outputAssocCollection{"MCRecoParticleAssociation", - Gaudi::DataHandle::Writer, this}; + const auto& mc = *mc_ptr; + const auto& tracks = *tracks_ptr; + auto& part = *part_ptr; + auto& assoc = *assoc_ptr; - // Matching momentum tolerance requires 10% by default; - Gaudi::Property m_pRelativeTolerance{this, "pRelativeTolerance", {0.1}}; - // Matching phi tolerance of 10 mrad - Gaudi::Property m_phiTolerance{this, "phiTolerance", {0.030}}; - // Matchin eta tolerance of 0.1 - Gaudi::Property m_etaTolerance{this, "etaTolerance", {0.2}}; - -public: - ParticlesWithTruthPID(const std::string& name, ISvcLocator* svcLoc) : GaudiAlgorithm(name, svcLoc) { - declareProperty("inputMCParticles", m_inputTruthCollection, "MCParticles"); - declareProperty("inputTrackParameters", m_inputTrackCollection, "outputTrackParameters"); - declareProperty("outputParticles", m_outputParticleCollection, "ReconstructedParticles"); - declareProperty("outputAssociations", m_outputAssocCollection, "MCRecoParticleAssociation"); - } - StatusCode initialize() override { - if (GaudiAlgorithm::initialize().isFailure()) { - return StatusCode::FAILURE; - } - return StatusCode::SUCCESS; - } - StatusCode execute() override { - // input collection - const auto& mc = *(m_inputTruthCollection.get()); - const auto& tracks = *(m_inputTrackCollection.get()); - auto& part = *(m_outputParticleCollection.createAndPut()); - auto& assoc = *(m_outputAssocCollection.createAndPut()); - - const double sinPhiOver2Tolerance = sin(0.5 * m_phiTolerance); - std::vector consumed(mc.size(), false); - for (const auto& trk : tracks) { - const auto mom = edm4eic::sphericalToVector(1.0 / std::abs(trk.getQOverP()), trk.getTheta(), trk.getPhi()); - const auto charge_rec = trk.getCharge(); - // utility variables for matching - int best_match = -1; - double best_delta = std::numeric_limits::max(); - for (size_t ip = 0; ip < mc.size(); ++ip) { - const auto& mcpart = mc[ip]; - if (consumed[ip] || mcpart.getGeneratorStatus() > 1 || mcpart.getCharge() == 0 || - mcpart.getCharge() * charge_rec < 0) { - if (msgLevel(MSG::DEBUG)) { - debug() << "ignoring non-primary/neutral/opposite charge particle" << endmsg; - } - continue; + const double sinPhiOver2Tolerance = sin(0.5 * m_phiTolerance); + std::vector consumed(mc.size(), false); + for (const auto& trk : tracks) { + const auto mom = + edm4eic::sphericalToVector(1.0 / std::abs(trk.getQOverP()), trk.getTheta(), trk.getPhi()); + const auto charge_rec = trk.getCharge(); + // utility variables for matching + int best_match = -1; + double best_delta = std::numeric_limits::max(); + for (size_t ip = 0; ip < mc.size(); ++ip) { + const auto& mcpart = mc[ip]; + if (consumed[ip] || mcpart.getGeneratorStatus() > 1 || mcpart.getCharge() == 0 || + mcpart.getCharge() * charge_rec < 0) { + if (aboveDebugThreshold()) { + debug() << "ignoring non-primary/neutral/opposite charge particle" << endmsg; } - const auto& p = mcpart.getMomentum(); - const auto p_mag = std::hypot(p.x, p.y, p.z); - const auto p_phi = std::atan2(p.y, p.x); - const auto p_eta = std::atanh(p.z / p_mag); - const double dp_rel = std::abs((edm4eic::magnitude(mom) - p_mag) / p_mag); - // check the tolerance for sin(dphi/2) to avoid the hemisphere problem and allow - // for phi rollovers - const double dsphi = std::abs(sin(0.5 * (edm4eic::angleAzimuthal(mom) - p_phi))); - const double deta = std::abs((edm4eic::eta(mom) - p_eta)); + continue; + } + const auto& p = mcpart.getMomentum(); + const auto p_mag = std::hypot(p.x, p.y, p.z); + const auto p_phi = std::atan2(p.y, p.x); + const auto p_eta = std::atanh(p.z / p_mag); + const double dp_rel = std::abs((edm4eic::magnitude(mom) - p_mag) / p_mag); + // check the tolerance for sin(dphi/2) to avoid the hemisphere problem and allow + // for phi rollovers + const double dsphi = std::abs(sin(0.5 * (edm4eic::angleAzimuthal(mom) - p_phi))); + const double deta = std::abs((edm4eic::eta(mom) - p_eta)); - if (dp_rel < m_pRelativeTolerance && deta < m_etaTolerance && dsphi < sinPhiOver2Tolerance) { - const double delta = - std::hypot(dp_rel / m_pRelativeTolerance, deta / m_etaTolerance, dsphi / sinPhiOver2Tolerance); - if (delta < best_delta) { - best_match = ip; - best_delta = delta; - } + if (dp_rel < m_pRelativeTolerance && deta < m_etaTolerance && dsphi < sinPhiOver2Tolerance) { + const double delta = std::hypot(dp_rel / m_pRelativeTolerance, deta / m_etaTolerance, + dsphi / sinPhiOver2Tolerance); + if (delta < best_delta) { + best_match = ip; + best_delta = delta; } } - auto rec_part = part.create(); - int32_t best_pid = 0; - auto referencePoint = rec_part.referencePoint(); - // float time = 0; - float mass = 0; - if (best_match >= 0) { - consumed[best_match] = true; - const auto& mcpart = mc[best_match]; - best_pid = mcpart.getPDG(); - referencePoint = { - static_cast(mcpart.getVertex().x), static_cast(mcpart.getVertex().y), - static_cast(mcpart.getVertex().z)}; // @TODO: not sure if vertex/reference poitn makes sense here - // time = mcpart.getTime(); - mass = mcpart.getMass(); - } - rec_part.setType(static_cast(best_match >= 0 ? 0 : -1)); // @TODO: determine type codes - rec_part.setEnergy(std::hypot(edm4eic::magnitude(mom), mass)); - rec_part.setMomentum(mom); - rec_part.setReferencePoint(referencePoint); - rec_part.setCharge(charge_rec); - rec_part.setMass(mass); - rec_part.setGoodnessOfPID(1); // perfect PID - rec_part.setPDG(best_pid); - // rec_part.covMatrix() // @TODO: covariance matrix on 4-momentum - // Also write MC <--> truth particle association if match was found - if (best_match >= 0) { - auto rec_assoc = assoc.create(); - rec_assoc.setRecID(rec_part.getObjectID().index); - rec_assoc.setSimID(mc[best_match].getObjectID().index); - rec_assoc.setWeight(1); - rec_assoc.setRec(rec_part); - //rec_assoc.setSim(mc[best_match]); - } - if (msgLevel(MSG::DEBUG)) { - if (best_match > 0) { - const auto& mcpart = mc[best_match]; - debug() << fmt::format("Matched track with MC particle {}\n", best_match) << endmsg; - debug() << fmt::format(" - Track: (mom: {}, theta: {}, phi: {}, charge: {})", edm4eic::magnitude(mom), - edm4eic::anglePolar(mom), edm4eic::angleAzimuthal(mom), charge_rec) - << endmsg; - const auto& p = mcpart.getMomentum(); - const auto p_mag = edm4eic::magnitude(p); - const auto p_phi = edm4eic::angleAzimuthal(p); - const auto p_theta = edm4eic::anglePolar(p); - debug() << fmt::format(" - MC particle: (mom: {}, theta: {}, phi: {}, charge: {}, type: {}", p_mag, p_theta, - p_phi, mcpart.getCharge(), mcpart.getPDG()) - << endmsg; - } else { - debug() << fmt::format("Did not find a good match for track \n") << endmsg; - debug() << fmt::format(" - Track: (mom: {}, theta: {}, phi: {}, charge: {})", edm4eic::magnitude(mom), - edm4eic::anglePolar(mom), edm4eic::angleAzimuthal(mom), charge_rec) - << endmsg; - } + } + auto rec_part = part.create(); + int32_t best_pid = 0; + auto referencePoint = rec_part.referencePoint(); + // float time = 0; + float mass = 0; + if (best_match >= 0) { + consumed[best_match] = true; + const auto& mcpart = mc[best_match]; + best_pid = mcpart.getPDG(); + referencePoint = { + static_cast(mcpart.getVertex().x), static_cast(mcpart.getVertex().y), + static_cast( + mcpart.getVertex().z)}; // @TODO: not sure if vertex/reference poitn makes sense here + // time = mcpart.getTime(); + mass = mcpart.getMass(); + } + rec_part.setType(static_cast(best_match >= 0 ? 0 : -1)); // @TODO: determine type codes + rec_part.setEnergy(std::hypot(edm4eic::magnitude(mom), mass)); + rec_part.setMomentum(mom); + rec_part.setReferencePoint(referencePoint); + rec_part.setCharge(charge_rec); + rec_part.setMass(mass); + rec_part.setGoodnessOfPID(1); // perfect PID + rec_part.setPDG(best_pid); + // rec_part.covMatrix() // @TODO: covariance matrix on 4-momentum + // Also write MC <--> truth particle association if match was found + if (best_match >= 0) { + auto rec_assoc = assoc.create(); + rec_assoc.setRecID(rec_part.getObjectID().index); + rec_assoc.setSimID(mc[best_match].getObjectID().index); + rec_assoc.setWeight(1); + rec_assoc.setRec(rec_part); + // rec_assoc.setSim(mc[best_match]); + } + if (aboveDebugThreshold()) { + if (best_match > 0) { + const auto& mcpart = mc[best_match]; + debug() << fmt::format("Matched track with MC particle {}\n", best_match) << endmsg; + debug() << fmt::format(" - Track: (mom: {}, theta: {}, phi: {}, charge: {})", + edm4eic::magnitude(mom), edm4eic::anglePolar(mom), + edm4eic::angleAzimuthal(mom), charge_rec) + << endmsg; + const auto& p = mcpart.getMomentum(); + const auto p_mag = edm4eic::magnitude(p); + const auto p_phi = edm4eic::angleAzimuthal(p); + const auto p_theta = edm4eic::anglePolar(p); + debug() << fmt::format( + " - MC particle: (mom: {}, theta: {}, phi: {}, charge: {}, type: {}", p_mag, + p_theta, p_phi, mcpart.getCharge(), mcpart.getPDG()) + << endmsg; + } else { + debug() << fmt::format("Did not find a good match for track \n") << endmsg; + debug() << fmt::format(" - Track: (mom: {}, theta: {}, phi: {}, charge: {})", + edm4eic::magnitude(mom), edm4eic::anglePolar(mom), + edm4eic::angleAzimuthal(mom), charge_rec) + << endmsg; } } - - return StatusCode::SUCCESS; } -}; // namespace Jug::Fast - -// NOLINTNEXTLINE(cppcoreguidelines-avoid-non-const-global-variables) -DECLARE_COMPONENT(ParticlesWithTruthPID) +} -} // namespace Jug::Fast +} // namespace algorithms::truth