From 597b5691cd9940114ddc03c23d4be71030934388 Mon Sep 17 00:00:00 2001 From: Sylvester Joosten Date: Tue, 27 Sep 2022 23:20:50 +0000 Subject: [PATCH 01/11] Initial implementation, needs testing --- .../core/include/algorithms/random.h | 112 ++++++++++++++++++ external/algorithms/core/src/dummy.cpp | 1 + external/algorithms/core/src/random.cpp | 30 +++++ 3 files changed, 143 insertions(+) create mode 100644 external/algorithms/core/include/algorithms/random.h create mode 100644 external/algorithms/core/src/random.cpp diff --git a/external/algorithms/core/include/algorithms/random.h b/external/algorithms/core/include/algorithms/random.h new file mode 100644 index 00000000..719dd885 --- /dev/null +++ b/external/algorithms/core/include/algorithms/random.h @@ -0,0 +1,112 @@ +#pragma once + +#include +#include +#include +#include +#include +#include +#include + +#include +#include + +namespace algorithms { + +// Cached link to the underlying generator 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; +}; + +// thread-safe random generator +class RandomSvc : public LoggedService { +public: + using value_type = CachedBitGenerator::value_type; + using GenFunc = CachedBitGenerator::GenFunc; + + void init(); + void init(const GenFunc& gen); + +public: + class Generator { + public: + Generator(const GenFunc& gen, const size_t cache_size) : m_gen{gen, cache_size} {} + + template Int uniform_int(const Int min, const Int max) { + 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) { + std::uniform_real_distribution d{min, max}; + std::lock_guard lock{m_mutex}; + return d(m_gen); + } + template Int poisson(const Int mean) { + std::poisson_distribution d{mean}; + std::lock_guard lock{m_mutex}; + return d(m_gen); + } + template Float exponential(const Float lambda) { + std::exponential_distribution d{lambda}; + std::lock_guard lock{m_mutex}; + return d(m_gen); + } + template Float gaussian(const Float mu, const Float sigma) { + std::normal_distribution d{mu, sigma}; + std::lock_guard lock{m_mutex}; + return d(m_gen); + } + + private: + CachedBitGenerator m_gen; + std::mutex m_mutex; + }; + Generator generator() { return {m_gen, m_cache_size}; } + +private: + GenFunc createBitGenerator(const size_t seed = 0); + + GenFunc m_gen{createBitGenerator()}; + 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 a264023d..5041865f 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 00000000..c1e8fa75 --- /dev/null +++ b/external/algorithms/core/src/random.cpp @@ -0,0 +1,30 @@ +#include + +namespace algorithms { + +void RandomSvc::init() { + if (m_seed.hasValue()) { + info() << "Custom random seed requested: " << m_seed << endmsg; + m_gen = createBitGenerator(m_seed); + } +} +void RandomSvc::init(const RandomSvc::GenFunc& 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; + } +} +RandomSvc::GenFunc RandomSvc::createBitGenerator(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; + }; +} +} // namespace algorithms -- GitLab From 825508d95ea040f3fbd04ac1b0f73533368b212c Mon Sep 17 00:00:00 2001 From: Sylvester Joosten Date: Wed, 28 Sep 2022 22:31:00 +0000 Subject: [PATCH 02/11] Ported a second algorithm --- JugFast/CMakeLists.txt | 2 + .../src/components/ParticlesWithTruthPID.cpp | 161 +---------- external/algorithms/CMakeLists.txt | 2 +- .../algorithms/calorimetry/ClusterRecoCoG.h | 1 - external/algorithms/truth/CMakeLists.txt | 1 + .../truth/ParticlesWithTruthPID.cpp | 163 ------------ .../algorithms/truth/ParticlesWithTruthPID.h | 39 +++ .../truth/src/ParticlesWithTruthPID.cpp | 249 ++++++++---------- 8 files changed, 152 insertions(+), 466 deletions(-) delete mode 100644 external/algorithms/truth/include/algorithms/truth/ParticlesWithTruthPID.cpp create mode 100644 external/algorithms/truth/include/algorithms/truth/ParticlesWithTruthPID.h diff --git a/JugFast/CMakeLists.txt b/JugFast/CMakeLists.txt index 91b9c9b9..3b937ef3 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/ParticlesWithTruthPID.cpp b/JugFast/src/components/ParticlesWithTruthPID.cpp index 40177ae4..face0d07 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 393a3a02..68401fb7 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 b43090ac..48f920e9 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 diff --git a/external/algorithms/truth/CMakeLists.txt b/external/algorithms/truth/CMakeLists.txt index cca38a85..4695f160 100644 --- a/external/algorithms/truth/CMakeLists.txt +++ b/external/algorithms/truth/CMakeLists.txt @@ -12,6 +12,7 @@ set(TARGETS ${TARGETS} ${LIBRARY} PARENT_SCOPE) # FIXME: adding one by one #file(GLOB SRC CONFIGURE_DEPENDS src/*.cpp) set(SRC + src/ParticlesWithTruthPID.cpp ) add_library(${LIBRARY} SHARED ${SRC}) 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 40177ae4..00000000 --- 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 00000000..934b1cb5 --- /dev/null +++ b/external/algorithms/truth/include/algorithms/truth/ParticlesWithTruthPID.h @@ -0,0 +1,39 @@ +// SPDX-License-Identifier: LGPL-3.0-or-later +// Copyright (C) 2022 Sylvester Joosten, Wouter Deconinck + +#include +#include + +// Event Model related classes +#include "edm4eic/MCRecoParticleAssociationCollection.h" +#include "edm4eic/ReconstructedParticleCollection.h" +#include "edm4eic/TrackParametersCollection.h" +#include "edm4hep/MCParticleCollection.h" + +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&); + +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/ParticlesWithTruthPID.cpp b/external/algorithms/truth/src/ParticlesWithTruthPID.cpp index 40177ae4..dedb636c 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 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 -- GitLab From d00747c25c67009eb8117b8e25317bcb983ab37c Mon Sep 17 00:00:00 2001 From: Sylvester Joosten Date: Wed, 28 Sep 2022 22:48:01 +0000 Subject: [PATCH 03/11] Add MC2SmearedParticles --- .../calorimetry/include/algorithms/calorimetry/ClusterRecoCoG.h | 2 -- .../truth/{MC2SmearedParticle.cpp => MC2SmearedParticle.h} | 0 2 files changed, 2 deletions(-) rename external/algorithms/truth/include/algorithms/truth/{MC2SmearedParticle.cpp => MC2SmearedParticle.h} (100%) diff --git a/external/algorithms/calorimetry/include/algorithms/calorimetry/ClusterRecoCoG.h b/external/algorithms/calorimetry/include/algorithms/calorimetry/ClusterRecoCoG.h index 48f920e9..1d788663 100644 --- a/external/algorithms/calorimetry/include/algorithms/calorimetry/ClusterRecoCoG.h +++ b/external/algorithms/calorimetry/include/algorithms/calorimetry/ClusterRecoCoG.h @@ -33,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 diff --git a/external/algorithms/truth/include/algorithms/truth/MC2SmearedParticle.cpp b/external/algorithms/truth/include/algorithms/truth/MC2SmearedParticle.h similarity index 100% rename from external/algorithms/truth/include/algorithms/truth/MC2SmearedParticle.cpp rename to external/algorithms/truth/include/algorithms/truth/MC2SmearedParticle.h -- GitLab From f37b783f9f1c3958d51d9303800a577fd8dca1a0 Mon Sep 17 00:00:00 2001 From: Sylvester Joosten Date: Wed, 28 Sep 2022 23:08:52 +0000 Subject: [PATCH 04/11] Made RNG work nice with const-correctness --- .../core/include/algorithms/detail/random.h | 48 +++++++ .../core/include/algorithms/random.h | 123 ++++++---------- external/algorithms/core/src/random.cpp | 8 +- external/algorithms/truth/CMakeLists.txt | 1 + .../algorithms/truth/MC2SmearedParticle.h | 100 +++---------- .../algorithms/truth/ParticlesWithTruthPID.h | 9 +- .../truth/src/MC2SmearedParticle.cpp | 135 +++++++----------- 7 files changed, 171 insertions(+), 253 deletions(-) create mode 100644 external/algorithms/core/include/algorithms/detail/random.h 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 00000000..cc43d754 --- /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 index 719dd885..cad4a7fd 100644 --- a/external/algorithms/core/include/algorithms/random.h +++ b/external/algorithms/core/include/algorithms/random.h @@ -1,107 +1,78 @@ #pragma once -#include #include #include -#include #include #include -#include +#include #include #include namespace algorithms { -// Cached link to the underlying generator allowing for multiple instances to be evaluated -// in parallel. Specs: -// - GenFunc is required to return N random numbers between 0 and +// 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() -// - 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 { +// - 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: - 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} {} + Generator(const RandomEngineCB& gen, const size_t cache_size) : m_gen{gen, cache_size} {} - value_type operator()() { - if (m_index >= m_cache.size()) { - refresh(); - } - return m_cache[m_index++]; + 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); } - - 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; + 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); } - GenFunc m_gen; - std::vector m_cache; - size_t m_index; +private: + mutable detail::CachedBitGenerator m_gen; + mutable std::mutex m_mutex; }; -// thread-safe random generator +// 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 = CachedBitGenerator::value_type; - using GenFunc = CachedBitGenerator::GenFunc; + using value_type = detail::CachedBitGenerator::value_type; void init(); - void init(const GenFunc& gen); - -public: - class Generator { - public: - Generator(const GenFunc& gen, const size_t cache_size) : m_gen{gen, cache_size} {} - - template Int uniform_int(const Int min, const Int max) { - 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) { - std::uniform_real_distribution d{min, max}; - std::lock_guard lock{m_mutex}; - return d(m_gen); - } - template Int poisson(const Int mean) { - std::poisson_distribution d{mean}; - std::lock_guard lock{m_mutex}; - return d(m_gen); - } - template Float exponential(const Float lambda) { - std::exponential_distribution d{lambda}; - std::lock_guard lock{m_mutex}; - return d(m_gen); - } - template Float gaussian(const Float mu, const Float sigma) { - std::normal_distribution d{mu, sigma}; - std::lock_guard lock{m_mutex}; - return d(m_gen); - } - - private: - CachedBitGenerator m_gen; - std::mutex m_mutex; - }; + void init(const RandomEngineCB& gen); Generator generator() { return {m_gen, m_cache_size}; } private: - GenFunc createBitGenerator(const size_t seed = 0); + RandomEngineCB createEngine(const size_t seed = 0); - GenFunc m_gen{createBitGenerator()}; + RandomEngineCB m_gen{createEngine()}; Property m_seed{this, "seed"}; Property m_cache_size{this, "cacheSize", 1024}; std::mutex m_mutex; diff --git a/external/algorithms/core/src/random.cpp b/external/algorithms/core/src/random.cpp index c1e8fa75..9956b367 100644 --- a/external/algorithms/core/src/random.cpp +++ b/external/algorithms/core/src/random.cpp @@ -5,10 +5,10 @@ namespace algorithms { void RandomSvc::init() { if (m_seed.hasValue()) { info() << "Custom random seed requested: " << m_seed << endmsg; - m_gen = createBitGenerator(m_seed); + m_gen = createEngine(m_seed); } } -void RandomSvc::init(const RandomSvc::GenFunc& gen) { +void RandomSvc::init(const RandomEngineCB& gen) { init(); info() << "Loading external generator function." << endmsg; m_gen = gen; @@ -17,12 +17,12 @@ void RandomSvc::init(const RandomSvc::GenFunc& gen) { << endmsg; } } -RandomSvc::GenFunc RandomSvc::createBitGenerator(const size_t seed) { +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::vector ret{size}; std::generate(ret.begin(), ret.end(), gen); return ret; }; diff --git a/external/algorithms/truth/CMakeLists.txt b/external/algorithms/truth/CMakeLists.txt index 4695f160..7c465f92 100644 --- a/external/algorithms/truth/CMakeLists.txt +++ b/external/algorithms/truth/CMakeLists.txt @@ -12,6 +12,7 @@ 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 ) diff --git a/external/algorithms/truth/include/algorithms/truth/MC2SmearedParticle.h b/external/algorithms/truth/include/algorithms/truth/MC2SmearedParticle.h index 6968e37f..b0bda623 100644 --- a/external/algorithms/truth/include/algorithms/truth/MC2SmearedParticle.h +++ b/external/algorithms/truth/include/algorithms/truth/MC2SmearedParticle.h @@ -1,99 +1,31 @@ // 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 +#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*/}; +using MC2SmearedParticleAlgorithm = Algorithm, + Output>; +class MC2SmearedParticle : public MC2SmearedParticleAlgorithm { 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; - } + MC2SmearedParticle(std::string_view name) + : MC2SmearedParticleAlgorithm{name, {"inputMCParticles"}, {"outputParticles"}} {} - 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; - } + void init(); + void process(const Input&, const Output&) const; - // 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; +private: + Generator m_rng = RandomSvc::instance().generator(); - 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; - } + Property m_smearing{this, "smearing", 0.01 /* 1 percent*/}; }; -// NOLINTNEXTLINE(cppcoreguidelines-avoid-non-const-global-variables) -DECLARE_COMPONENT(MC2SmearedParticle) +} // namespace algorithms::truth -} // namespace Jug::Fast diff --git a/external/algorithms/truth/include/algorithms/truth/ParticlesWithTruthPID.h b/external/algorithms/truth/include/algorithms/truth/ParticlesWithTruthPID.h index 934b1cb5..c39894ee 100644 --- a/external/algorithms/truth/include/algorithms/truth/ParticlesWithTruthPID.h +++ b/external/algorithms/truth/include/algorithms/truth/ParticlesWithTruthPID.h @@ -2,13 +2,12 @@ // Copyright (C) 2022 Sylvester Joosten, Wouter Deconinck #include -#include // Event Model related classes -#include "edm4eic/MCRecoParticleAssociationCollection.h" -#include "edm4eic/ReconstructedParticleCollection.h" -#include "edm4eic/TrackParametersCollection.h" -#include "edm4hep/MCParticleCollection.h" +#include +#include +#include +#include namespace algorithms::truth { diff --git a/external/algorithms/truth/src/MC2SmearedParticle.cpp b/external/algorithms/truth/src/MC2SmearedParticle.cpp index 6968e37f..89d9ccd3 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 -- GitLab From 03598406c10a23b4529f9e1ce2fb0731b34b566a Mon Sep 17 00:00:00 2001 From: Sylvester Joosten Date: Wed, 28 Sep 2022 23:14:02 +0000 Subject: [PATCH 05/11] Update Juggler plugin to use algorithms --- JugFast/src/components/MC2SmearedParticle.cpp | 97 +------------------ 1 file changed, 3 insertions(+), 94 deletions(-) diff --git a/JugFast/src/components/MC2SmearedParticle.cpp b/JugFast/src/components/MC2SmearedParticle.cpp index 6968e37f..148c870c 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) -- GitLab From dd7e68134aa790ab42235165a4d41e5ba9220bd1 Mon Sep 17 00:00:00 2001 From: Sylvester Joosten Date: Wed, 28 Sep 2022 23:28:19 +0000 Subject: [PATCH 06/11] remove random stuff from cpp to h for now as the CMake setup isn't working properly --- .../core/include/algorithms/random.h | 36 +++++++++++++++++-- external/algorithms/core/src/random.cpp | 3 ++ 2 files changed, 37 insertions(+), 2 deletions(-) diff --git a/external/algorithms/core/include/algorithms/random.h b/external/algorithms/core/include/algorithms/random.h index cad4a7fd..b9b37236 100644 --- a/external/algorithms/core/include/algorithms/random.h +++ b/external/algorithms/core/include/algorithms/random.h @@ -65,13 +65,45 @@ 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); - Generator generator() { return {m_gen, m_cache_size}; } +#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 = 0); + 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}; diff --git a/external/algorithms/core/src/random.cpp b/external/algorithms/core/src/random.cpp index 9956b367..711e35ab 100644 --- a/external/algorithms/core/src/random.cpp +++ b/external/algorithms/core/src/random.cpp @@ -2,6 +2,8 @@ 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; @@ -27,4 +29,5 @@ RandomEngineCB RandomSvc::createEngine(const size_t seed) { return ret; }; } +#endif } // namespace algorithms -- GitLab From 95a77cb409a5b3fbf2f5a0f918fc28be1f1f2efb Mon Sep 17 00:00:00 2001 From: Sylvester Joosten Date: Wed, 28 Sep 2022 23:28:39 +0000 Subject: [PATCH 07/11] fix const-correctness --- .../include/algorithms/calorimetry/ClusterRecoCoG.h | 2 +- external/algorithms/calorimetry/src/ClusterRecoCoG.cpp | 2 +- .../truth/include/algorithms/truth/ParticlesWithTruthPID.h | 6 +++--- external/algorithms/truth/src/ParticlesWithTruthPID.cpp | 2 +- 4 files changed, 6 insertions(+), 6 deletions(-) diff --git a/external/algorithms/calorimetry/include/algorithms/calorimetry/ClusterRecoCoG.h b/external/algorithms/calorimetry/include/algorithms/calorimetry/ClusterRecoCoG.h index 1d788663..96502e07 100644 --- a/external/algorithms/calorimetry/include/algorithms/calorimetry/ClusterRecoCoG.h +++ b/external/algorithms/calorimetry/include/algorithms/calorimetry/ClusterRecoCoG.h @@ -42,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 95018716..a456cd3f 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/truth/include/algorithms/truth/ParticlesWithTruthPID.h b/external/algorithms/truth/include/algorithms/truth/ParticlesWithTruthPID.h index c39894ee..d00fb76c 100644 --- a/external/algorithms/truth/include/algorithms/truth/ParticlesWithTruthPID.h +++ b/external/algorithms/truth/include/algorithms/truth/ParticlesWithTruthPID.h @@ -19,11 +19,11 @@ class ParticlesWithTruthPID : public ParticlesWithTruthPIDAlgorithm { public: ParticlesWithTruthPID(std::string_view name) : ParticlesWithTruthPIDAlgorithm{name, - {"inputMCParticles", "inputTrackParameters"}, - {"outputParticles", "outputAssociations"}} {} + {"inputMCParticles", "inputTrackParameters"}, + {"outputParticles", "outputAssociations"}} {} void init(); - void process(const Input&, const Output&); + void process(const Input&, const Output&) const; private: // Matching momentum tolerance requires 10% by default; diff --git a/external/algorithms/truth/src/ParticlesWithTruthPID.cpp b/external/algorithms/truth/src/ParticlesWithTruthPID.cpp index dedb636c..4cbe18e2 100644 --- a/external/algorithms/truth/src/ParticlesWithTruthPID.cpp +++ b/external/algorithms/truth/src/ParticlesWithTruthPID.cpp @@ -15,7 +15,7 @@ void ParticlesWithTruthPID::init() { ; // do nothing } void ParticlesWithTruthPID::process(const ParticlesWithTruthPID::Input& input, - const ParticlesWithTruthPID::Output& output) { + const ParticlesWithTruthPID::Output& output) const { const auto [mc_ptr, tracks_ptr] = input; auto [part_ptr, assoc_ptr] = output; -- GitLab From 56e19fea96955f86b8f36add3405c3a9ce07a8bd Mon Sep 17 00:00:00 2001 From: Sylvester Joosten Date: Wed, 28 Sep 2022 23:40:31 +0000 Subject: [PATCH 08/11] Make AlgoServiceSvc aware of algorithms::RandomSvc (but don't actually setup for now) --- JugAlgo/src/components/AlgoServiceSvc.cpp | 5 +++++ 1 file changed, 5 insertions(+) diff --git a/JugAlgo/src/components/AlgoServiceSvc.cpp b/JugAlgo/src/components/AlgoServiceSvc.cpp index 4af79335..36a90b80 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 @@ -83,6 +84,10 @@ 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() << "No setup for algorithms::RandomSvc --> using internal STL 64-bit MT engine" + << endmsg; } else { fatal() << "Unknown service encountered, please implement the necessary framework hooks" << endmsg; -- GitLab From 61045c5ea3671d349b26f4e52a3b1915f7d9b520 Mon Sep 17 00:00:00 2001 From: Sylvester Joosten Date: Thu, 29 Sep 2022 14:43:08 +0000 Subject: [PATCH 09/11] fix argument name --- .../truth/include/algorithms/truth/MC2SmearedParticle.h | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/external/algorithms/truth/include/algorithms/truth/MC2SmearedParticle.h b/external/algorithms/truth/include/algorithms/truth/MC2SmearedParticle.h index b0bda623..e7ee0dfb 100644 --- a/external/algorithms/truth/include/algorithms/truth/MC2SmearedParticle.h +++ b/external/algorithms/truth/include/algorithms/truth/MC2SmearedParticle.h @@ -16,7 +16,7 @@ using MC2SmearedParticleAlgorithm = Algorithm Date: Thu, 29 Sep 2022 19:04:38 +0000 Subject: [PATCH 10/11] Use propper constructor for a std::vector of size N --- external/algorithms/core/include/algorithms/random.h | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/external/algorithms/core/include/algorithms/random.h b/external/algorithms/core/include/algorithms/random.h index b9b37236..6826aada 100644 --- a/external/algorithms/core/include/algorithms/random.h +++ b/external/algorithms/core/include/algorithms/random.h @@ -96,7 +96,7 @@ private: static std::mutex m; static std::mt19937_64 gen{seed}; std::lock_guard lock{m}; - std::vector ret{size}; + std::vector ret(size); std::generate(ret.begin(), ret.end(), gen); return ret; }; -- GitLab From 0daad6488911b67d35299e87f94d60255fbffec3 Mon Sep 17 00:00:00 2001 From: Sylvester Joosten Date: Thu, 29 Sep 2022 19:42:21 +0000 Subject: [PATCH 11/11] Make random seed configurable from the Gaudi end --- JugAlgo/src/components/AlgoServiceSvc.cpp | 9 +++++++-- 1 file changed, 7 insertions(+), 2 deletions(-) diff --git a/JugAlgo/src/components/AlgoServiceSvc.cpp b/JugAlgo/src/components/AlgoServiceSvc.cpp index 36a90b80..90cb775c 100644 --- a/JugAlgo/src/components/AlgoServiceSvc.cpp +++ b/JugAlgo/src/components/AlgoServiceSvc.cpp @@ -23,6 +23,7 @@ public: private: SmartIF m_geoSvc; + Gaudi::Property m_randomSeed{this, "randomSeed", 1}; }; DECLARE_COMPONENT(AlgoServiceSvc) @@ -86,8 +87,12 @@ StatusCode AlgoServiceSvc::initialize() { geo->init(m_geoSvc->detector()); } else if (name == algorithms::RandomSvc::kName) { // setup random service - info() << "No setup for algorithms::RandomSvc --> using internal STL 64-bit MT engine" - << endmsg; + 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; -- GitLab