diff --git a/.clang-format b/.clang-format index 5737aca1ec7df306124ffc6d805befe3ff6ade12..5c443be69c45337924324a1f56fad9721dac62f6 100644 --- a/.clang-format +++ b/.clang-format @@ -9,7 +9,7 @@ SpaceAfterControlStatementKeyword: true PointerBindsToType: true IncludeBlocks: Preserve UseTab: Never -ColumnLimit: 120 +ColumnLimit: 100 NamespaceIndentation: Inner AlignConsecutiveAssignments: true ... diff --git a/CMakeLists.txt b/CMakeLists.txt index d385e94f9b1f6a62ae5d31ac6657399611ca2fc9..d8d743f7ff65307c10d1849c2413193aa1ea2f80 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -33,7 +33,7 @@ endif() # Set default build type set(default_build_type "Release") if(EXISTS "${CMAKE_SOURCE_DIR}/.git") - set(default_build_type "Debug") + set(default_build_type "RelWithDebInfo") endif() if(NOT CMAKE_BUILD_TYPE AND NOT CMAKE_CONFIGURATION_TYPES) message(STATUS "Setting build type to '${default_build_type}' as none was specified.") @@ -73,8 +73,12 @@ add_definitions("-DActs_VERSION_MAJOR=${Acts_VERSION_MAJOR}") add_definitions("-DActs_VERSION_MINOR=${Acts_VERSION_MINOR}") add_definitions("-DActs_VERSION_PATCH=${Acts_VERSION_PATCH}") +## algorithms dependency +add_subdirectory(external/algorithms) + find_package(Gaudi) add_subdirectory(JugBase) +add_subdirectory(JugAlgo) add_subdirectory(JugDigi) add_subdirectory(JugFast) add_subdirectory(JugPID) @@ -82,6 +86,7 @@ add_subdirectory(JugReco) add_subdirectory(JugTrack) gaudi_install(CMAKE) + # create and install Juggler.xenv file as it still has a use-case # TODO: update workflow to not need xenv files anymore include(cmake/xenv.cmake) diff --git a/JugAlgo/CMakeLists.txt b/JugAlgo/CMakeLists.txt new file mode 100644 index 0000000000000000000000000000000000000000..d97bace08bf6f5bb0ac4d6b93a29536362f2f2bb --- /dev/null +++ b/JugAlgo/CMakeLists.txt @@ -0,0 +1,51 @@ +# SPDX-License-Identifier: LGPL-3.0-or-later +# Copyright (C) 2022 Sylvester Joosten + +################################################################################ +# Package: JugAlgo +################################################################################ + +file(GLOB JugAlgo_sources CONFIGURE_DEPENDS src/*.cpp) +gaudi_add_library(JugAlgo + SOURCES + ${JugAlgo_sources} + LINK + Gaudi::GaudiKernel Gaudi::GaudiAlgLib + algocore + JugBase +) + +target_include_directories(JugAlgo PUBLIC + $ + $ +) + +target_compile_options(JugAlgo PRIVATE -Wno-suggest-override) + +file(GLOB JugAlgoPlugins_sources CONFIGURE_DEPENDS src/components/*.cpp) +gaudi_add_module(JugAlgoPlugins + SOURCES + ${JugAlgoPlugins_sources} + LINK + Gaudi::GaudiKernel Gaudi::GaudiAlgLib + JugBase + JugAlgo + algocore +) + +target_include_directories(JugAlgoPlugins PUBLIC + $ + $ +) + +target_compile_options(JugAlgoPlugins PRIVATE -Wno-suggest-override) + +install(TARGETS JugAlgo JugAlgoPlugins + EXPORT JugAlgoTargets + RUNTIME DESTINATION "${CMAKE_INSTALL_BINDIR}" COMPONENT bin + LIBRARY DESTINATION "${CMAKE_INSTALL_LIBDIR}" COMPONENT shlib + COMPONENT dev) + +if(BUILD_TESTING) + enable_testing() +endif() diff --git a/JugAlgo/JugAlgo/Algorithm.h b/JugAlgo/JugAlgo/Algorithm.h new file mode 100644 index 0000000000000000000000000000000000000000..937b3725962e2178afaffbb45271a52181c7a07f --- /dev/null +++ b/JugAlgo/JugAlgo/Algorithm.h @@ -0,0 +1,94 @@ +// SPDX-License-Identifier: LGPL-3.0-or-later +// Copyright (C) 2022 Sylvester Joosten + +#include +#include + +#include +#include + +#include +#include +#include +#include + +namespace Jug::Algo { + +template class Algorithm : public GaudiAlgorithm { +public: + using algo_type = AlgoImpl; + using input_type = typename algo_type::input_type; + using output_type = typename algo_type::output_type; + using Input = typename algo_type::Input; + using Output = typename algo_type::Output; + + Algorithm(const std::string& name, ISvcLocator* svcLoc) + : GaudiAlgorithm(name, svcLoc) + , m_algo{name} + , m_output{this, m_algo.outputNames()} + , m_input{this, m_algo.inputNames()} {} + + StatusCode initialize() override { + debug() << "Initializing " << name() << endmsg; + + // Grab the AlgoServiceSvc + m_algo_svc = service("AlgoServiceSvc"); + if (!m_algo_svc) { + error() << "Unable to get an instance of the AlgoServiceSvc" << endmsg; + return StatusCode::FAILURE; + } + + // Forward the log level of this algorithm + const algorithms::LogLevel level{ + static_cast(msgLevel() > 0 ? msgLevel() - 1 : 0)}; + debug() << "Setting the logger level to " << algorithms::logLevelName(level) << endmsg; + m_algo.level(level); + + // Init our data structures + debug() << "Initializing data structures" << endmsg; + m_input.init(); + m_output.init(); + + // call configure function that passes properties + debug() << "Configuring properties" << endmsg; + auto sc = configure(); + if (sc != StatusCode::SUCCESS) { + return sc; + } + + // call the internal algorithm init + debug() << "Initializing underlying algorithm " << m_algo.name() << endmsg; + m_algo.init(); + return StatusCode::SUCCESS; + } + + StatusCode execute() override { + try { + m_algo.process(m_input.get(), m_output.get()); + } catch (const std::exception& e) { + error() << e.what() << endmsg; + return StatusCode::FAILURE; + } + return StatusCode::SUCCESS; + } + + virtual StatusCode configure() = 0; + +protected: + template void setAlgoProp(std::string_view name, T&& value) { + m_algo.template setProperty(name, value); + } + template T getAlgoProp(std::string name) const { + return m_algo.template getProperty(name); + } + bool hasAlgoProp(std::string_view name) const { return m_algo.hasProperty(name); } + +private: + algo_type m_algo; + SmartIF m_algo_svc; + detail::DataProxy m_output; + detail::DataProxy m_input; +}; + +} // namespace Jug::Algo + diff --git a/JugAlgo/JugAlgo/IAlgoServiceSvc.h b/JugAlgo/JugAlgo/IAlgoServiceSvc.h new file mode 100644 index 0000000000000000000000000000000000000000..bc0c89b51b5bf5c8fe7345064770de04b85b9617 --- /dev/null +++ b/JugAlgo/JugAlgo/IAlgoServiceSvc.h @@ -0,0 +1,16 @@ +// SPDX-License-Identifier: LGPL-3.0-or-later +// Copyright (C) 2022 Sylvester Joosten + +#pragma once + +#include + +// Juggler bindings for all required algorithms services +// Will setup all necessary services (using the algorithms::ServiceSvc) + +class GAUDI_API IAlgoServiceSvc : virtual public IService { +public: + DeclareInterfaceID(IAlgoServiceSvc, 1, 0); + virtual ~IAlgoServiceSvc() {} + // No actual API needed, as this service will do all the magic behind the screens +}; diff --git a/JugAlgo/JugAlgo/detail/DataProxy.h b/JugAlgo/JugAlgo/detail/DataProxy.h new file mode 100644 index 0000000000000000000000000000000000000000..77b77c345aee13ea4b94c12487a04505839967e2 --- /dev/null +++ b/JugAlgo/JugAlgo/detail/DataProxy.h @@ -0,0 +1,156 @@ +#pragma once + +#include +#include +#include + +#include +#include + +#include +#include + +namespace Jug::Algo::detail { + +enum class DataMode : unsigned { kInput, kOutput }; + +// Generate properties for each of the data arguments +template class DataElement { + +public: + using value_type = std::conditional_t, + algorithms::output_type_t>; + using data_type = algorithms::data_type_t; + constexpr static const bool kIsOptional = algorithms::is_optional_v; + + template + DataElement(Owner* owner, std::string_view name) + : m_data_name{std::make_unique>(owner, std::string(name), "")} + , m_owner{owner} {} + void init() { + if (m_handle) { + // treat error: already initialized + } + if (!m_data_name->empty()) { + m_handle = std::make_unique>( + *m_data_name, + ((kMode == DataMode::kInput) ? Gaudi::DataHandle::Reader : Gaudi::DataHandle::Writer), + m_owner); + } else if (!algorithms::is_optional_v) { + // treat error: member not optional but no collection name given + } + } + value_type get() const { + if constexpr (kIsOptional) { + if (!m_handle) { + return nullptr; + } + } + if constexpr (kMode == DataMode::kInput) { + return m_handle->get(); + } else { + return m_handle->createAndPut(); + } + } + +private: + std::unique_ptr> + m_data_name; // This needs to be a pointer, else things go wrong once we go through + // createElements - probably something about passing the Property through an + // rvalue (or copy) constructor + std::unique_ptr> m_handle; + gsl::not_null m_owner; +}; + +// Specialization for vectors +template class DataElement, kMode> { +public: + using value_type = std::conditional_t, + algorithms::output_type_t>; + using data_type = algorithms::data_type_t; + + template + DataElement(Owner* owner, std::string_view name) + : m_data_names{std::make_unique>>( + owner, std::string(name), {})} + , m_owner{owner} {} + void init() { + if (!m_handles.empty()) { + // treat error: already initialized + } + if (!m_data_names->empty()) { + for (const auto& name : *m_data_names) { + if (!name.empty()) { + m_handles.emplace_back(std::make_unique>( + name, + (kMode == DataMode::kInput ? Gaudi::DataHandle::Reader : Gaudi::DataHandle::Writer), + m_owner)); + } else { + // treat error: empty name + } + } + } else { + // OK if nothing given here, no error + } + } + std::vector get() const { + std::vector ret; + for (auto& handle : m_handles) { + if constexpr (kMode == DataMode::kInput) { + ret.emplace_back(handle->get()); + } else { + ret.emplace_back(handle->createAndPut()); + } + } + return ret; + } + +private: + std::unique_ptr>> m_data_names; + std::vector>> m_handles; + gsl::not_null m_owner; +}; + +template +auto createElements(Owner* owner, const NamesArray& names, const Tuple&, std::index_sequence) + -> std::tuple, kMode>...> { + return {DataElement, kMode>(owner, std::get(names))...}; +} + +// Call ::get() on each element of the HandleTuple, and return the result in the format of +// ReturnTuple +template +ReturnTuple getElements(HandleTuple& handles, std::index_sequence) { + return {std::get(handles).get()...}; +} + +// Create data handle structure for all members + +template class DataProxy { +public: + static constexpr DataMode kMode = + (algorithms::is_input_v ? DataMode::kInput : DataMode::kOutput); + using value_type = typename Data::value_type; + using data_type = typename Data::data_type; + constexpr static size_t kSize = Data::kSize; + using names_type = typename Data::key_type; + using elements_type = + decltype(createElements(std::declval(), names_type(), data_type(), + std::make_index_sequence())); + + template + DataProxy(Owner* owner, const names_type& names) + : m_elements{ + createElements(owner, names, data_type(), std::make_index_sequence())} {} + void init() { + std::apply([](auto&&... el) { (el.init(), ...); }, m_elements); + } + value_type get() const { + return getElements(m_elements, std::make_index_sequence()); + } + +private: + elements_type m_elements; +}; + +} // namespace Jug::Algo::detail diff --git a/JugAlgo/README.md b/JugAlgo/README.md new file mode 100644 index 0000000000000000000000000000000000000000..2a48f8f603bac2902907e1264eb18e16411a62f7 --- /dev/null +++ b/JugAlgo/README.md @@ -0,0 +1 @@ +Juggler bindings for the `algorithms` library. diff --git a/JugAlgo/src/components/AlgoServiceSvc.cpp b/JugAlgo/src/components/AlgoServiceSvc.cpp new file mode 100644 index 0000000000000000000000000000000000000000..47bc6f49c819a79dab01d65519422f100227515a --- /dev/null +++ b/JugAlgo/src/components/AlgoServiceSvc.cpp @@ -0,0 +1,87 @@ +// SPDX-License-Identifier: LGPL-3.0-or-later +// Copyright (C) 2022 Sylvester Joosten + +#include + +#include +#include +#include + +#include +#include +#include + +// too many services? :P +class AlgoServiceSvc : public extends { +public: + AlgoServiceSvc(const std::string& name, ISvcLocator* svc) : base_class(name, svc) {} + virtual ~AlgoServiceSvc() = default; + + virtual StatusCode initialize() final; + virtual StatusCode finalize() final { return StatusCode::SUCCESS; } + +private: + SmartIF m_geoSvc; +}; + +DECLARE_COMPONENT(AlgoServiceSvc) + +// Implementation + +StatusCode AlgoServiceSvc::initialize() { + StatusCode sc = Service::initialize(); + if (!sc.isSuccess()) { + fatal() << "Error initializing AlgoServiceSvc" << endmsg; + return sc; + } + + auto& serviceSvc = algorithms::ServiceSvc::instance(); + info() << "ServiceSvc declared " << serviceSvc.services().size() << " services" << endmsg; + // loop over all services and handle each properly + // Note: this code is kind of dangerous, as getting the types wrong will lead to + // undefined runtime behavior. + for (auto [name, svc] : serviceSvc.services()) { + if (name == algorithms::LogSvc::kName) { + auto* logger = static_cast(svc); + const algorithms::LogLevel level{ + static_cast(msgLevel() > 0 ? msgLevel() - 1 : 0)}; + info() << "Setting up algorithms::LogSvc with default level " + << algorithms::logLevelName(level) << endmsg; + logger->defaultLevel(level); + logger->action( + [this](const algorithms::LogLevel l, std::string_view caller, std::string_view msg) { + const std::string text = fmt::format("[{}] {}", caller, msg); + if (l == algorithms::LogLevel::kCritical) { + this->fatal() << text << endmsg; + } else if (l == algorithms::LogLevel::kError) { + this->error() << text << endmsg; + } else if (l == algorithms::LogLevel::kWarning) { + this->warning() << text << endmsg; + } else if (l == algorithms::LogLevel::kInfo) { + this->info() << text << endmsg; + } else if (l == algorithms::LogLevel::kDebug) { + this->debug() << text << endmsg; + } else if (l == algorithms::LogLevel::kTrace) { + this->verbose() << text << endmsg; + } + }); + // set own log level to verbose so we actually display everything that is requested + // (this overrides what was initally set through the OutputLevel property) + updateMsgStreamOutputLevel(MSG::VERBOSE); + } else if (name == algorithms::GeoSvc::kName) { + // Setup geometry service + m_geoSvc = service("GeoSvc"); + if (!m_geoSvc) { + error() << "Unable to locate Geometry Service. " + << "Make sure you have GeoSvc in the right order in the configuration." << endmsg; + return StatusCode::FAILURE; + } + info() << "Setting up algorithms::GeoSvc" << endmsg; + auto* geo = static_cast(svc); + geo->init(m_geoSvc->detector()); + } + } + + info() << "AlgoServiceSvc initialized successfully" << endmsg; + return StatusCode::SUCCESS; +} diff --git a/JugAlgo/src/dummy.cpp b/JugAlgo/src/dummy.cpp new file mode 100644 index 0000000000000000000000000000000000000000..959f35252f94806f9ae198528c2e78c20722352a --- /dev/null +++ b/JugAlgo/src/dummy.cpp @@ -0,0 +1,6 @@ +#include +//#include + +namespace { +constexpr int doNothing() { return 1; } +} // namespace diff --git a/JugDigi/CMakeLists.txt b/JugDigi/CMakeLists.txt index 41d082b7ac5d9cb04e803f478fbccf499ce3cdf4..a177c3973eb3c5178422af7b4dd6349eaa8bb49e 100644 --- a/JugDigi/CMakeLists.txt +++ b/JugDigi/CMakeLists.txt @@ -11,7 +11,8 @@ gaudi_add_module(JugDigiPlugins ${JugDigiPlugins_sources} LINK Gaudi::GaudiAlgLib Gaudi::GaudiKernel - JugBase + JugBase JugAlgo + algocore algocalorimetry ROOT::Core ROOT::RIO ROOT::Tree EDM4HEP::edm4hep EDM4EIC::edm4eic diff --git a/JugDigi/src/components/CalorimeterHitDigi.cpp b/JugDigi/src/components/CalorimeterHitDigi.cpp index a108d7656b022ab4fad08521a3c12895a00a2e29..c93bbf2b4aab9ed6c8f8a349b747c03fe10037b9 100644 --- a/JugDigi/src/components/CalorimeterHitDigi.cpp +++ b/JugDigi/src/components/CalorimeterHitDigi.cpp @@ -1,264 +1,65 @@ // SPDX-License-Identifier: LGPL-3.0-or-later // Copyright (C) 2022 Chao Peng, Wouter Deconinck, Sylvester Joosten -// A general digitization for CalorimeterHit from simulation -// 1. Smear energy deposit with a/sqrt(E/GeV) + b + c/E or a/sqrt(E/GeV) (relative value) -// 2. Digitize the energy with dynamic ADC range and add pedestal (mean +- sigma) -// 3. Time conversion with smearing resolution (absolute value) -// 4. Signal is summed if the SumFields are provided -// -// Author: Chao Peng -// Date: 06/02/2021 +#include +#include -#include -#include -#include - -#include "GaudiAlg/GaudiTool.h" -#include "GaudiAlg/Transformer.h" -#include "GaudiKernel/PhysicalConstants.h" #include "Gaudi/Property.h" -#include "GaudiKernel/RndmGenerators.h" - -#include "DDRec/CellIDPositionConverter.h" -#include "DDSegmentation/BitFieldCoder.h" - -#include "JugBase/IGeoSvc.h" -#include "JugBase/DataHandle.h" - -#include "fmt/format.h" -#include "fmt/ranges.h" - -// Event Model related classes -#include "edm4hep/SimCalorimeterHitCollection.h" -#include "edm4eic/RawCalorimeterHitCollection.h" -#include "edm4eic/RawCalorimeterHitData.h" - - -using namespace Gaudi::Units; namespace Jug::Digi { - /** Generic calorimeter hit digitiziation. - * - * \ingroup digi - * \ingroup calorimetry - */ - class CalorimeterHitDigi : public GaudiAlgorithm { - public: - // additional smearing resolutions - Gaudi::Property> u_eRes{this, "energyResolutions", {}}; // a/sqrt(E/GeV) + b + c/(E/GeV) - Gaudi::Property m_tRes{this, "timeResolution", 0.0 * ns}; - // single hit energy deposition threshold - Gaudi::Property m_threshold{this, "threshold", 1. * keV}; - - // digitization settings - Gaudi::Property m_capADC{this, "capacityADC", 8096}; - Gaudi::Property m_dyRangeADC{this, "dynamicRangeADC", 100 * MeV}; - Gaudi::Property m_pedMeanADC{this, "pedestalMean", 400}; - Gaudi::Property m_pedSigmaADC{this, "pedestalSigma", 3.2}; - Gaudi::Property m_resolutionTDC{this, "resolutionTDC", 10 * ps}; - - Gaudi::Property m_corrMeanScale{this, "scaleResponse", 1.0}; - // These are better variable names for the "energyResolutions" array which is a bit - // magic @FIXME - //Gaudi::Property m_corrSigmaCoeffE{this, "responseCorrectionSigmaCoeffE", 0.0}; - //Gaudi::Property m_corrSigmaCoeffSqrtE{this, "responseCorrectionSigmaCoeffSqrtE", 0.0}; - - // signal sums - // @TODO: implement signal sums with timing - // field names to generate id mask, the hits will be grouped by masking the field - Gaudi::Property> u_fields{this, "signalSumFields", {}}; - // ref field ids are used for the merged hits, 0 is used if nothing provided - Gaudi::Property> u_refs{this, "fieldRefNumbers", {}}; - Gaudi::Property m_geoSvcName{this, "geoServiceName", "GeoSvc"}; - Gaudi::Property m_readout{this, "readoutClass", ""}; - - // unitless counterparts of inputs - double dyRangeADC{0}, stepTDC{0}, tRes{0}, eRes[3] = {0., 0., 0.}; - Rndm::Numbers m_normDist; - SmartIF m_geoSvc; - uint64_t id_mask{0}, ref_mask{0}; - - DataHandle m_inputHitCollection{ - "inputHitCollection", Gaudi::DataHandle::Reader, this}; - DataHandle m_outputHitCollection{ - "outputHitCollection", Gaudi::DataHandle::Writer, this}; - - // ill-formed: using GaudiAlgorithm::GaudiAlgorithm; - CalorimeterHitDigi(const std::string& name, ISvcLocator* svcLoc) - : GaudiAlgorithm(name, svcLoc) { - declareProperty("inputHitCollection", m_inputHitCollection, ""); - declareProperty("outputHitCollection", m_outputHitCollection, ""); - } - - StatusCode initialize() override - { - if (GaudiAlgorithm::initialize().isFailure()) { - return StatusCode::FAILURE; - } - // random number generator from service - auto randSvc = svc("RndmGenSvc", true); - auto sc = m_normDist.initialize(randSvc, Rndm::Gauss(0.0, 1.0)); - if (!sc.isSuccess()) { - return StatusCode::FAILURE; - } - // set energy resolution numbers - for (size_t i = 0; i < u_eRes.size() && i < 3; ++i) { - eRes[i] = u_eRes[i]; - } - - // using juggler internal units (GeV, mm, radian, ns) - dyRangeADC = m_dyRangeADC.value() / GeV; - tRes = m_tRes.value() / ns; - stepTDC = ns / m_resolutionTDC.value(); - - // need signal sum - if (!u_fields.value().empty()) { - m_geoSvc = service(m_geoSvcName); - // sanity checks - if (!m_geoSvc) { - error() << "Unable to locate Geometry Service. " - << "Make sure you have GeoSvc and SimSvc in the right order in the configuration." - << endmsg; - return StatusCode::FAILURE; - } - if (m_readout.value().empty()) { - error() << "readoutClass is not provided, it is needed to know the fields in readout ids" - << endmsg; - return StatusCode::FAILURE; - } - - // get decoders - try { - auto id_desc = m_geoSvc->detector()->readout(m_readout).idSpec(); - id_mask = 0; - std::vector> ref_fields; - for (size_t i = 0; i < u_fields.size(); ++i) { - id_mask |= id_desc.field(u_fields[i])->mask(); - // use the provided id number to find ref cell, or use 0 - int ref = i < u_refs.size() ? u_refs[i] : 0; - ref_fields.emplace_back(u_fields[i], ref); - } - ref_mask = id_desc.encode(ref_fields); - // debug() << fmt::format("Referece id mask for the fields {:#064b}", ref_mask) << endmsg; - } catch (...) { - error() << "Failed to load ID decoder for " << m_readout << endmsg; - return StatusCode::FAILURE; - } - id_mask = ~id_mask; - info() << fmt::format("ID mask in {:s}: {:#064b}", m_readout, id_mask) << endmsg; - return StatusCode::SUCCESS; - } - - return StatusCode::SUCCESS; - } - - StatusCode execute() override - { - if (!u_fields.value().empty()) { - signal_sum_digi(); - } else { - single_hits_digi(); - } - return StatusCode::SUCCESS; - } - - private: - void single_hits_digi() { - // input collections - const auto* const simhits = m_inputHitCollection.get(); - // Create output collections - auto* rawhits = m_outputHitCollection.createAndPut(); - for (const auto& ahit : *simhits) { - // Note: juggler internal unit of energy is GeV - const double eDep = ahit.getEnergy(); - - // apply additional calorimeter noise to corrected energy deposit - const double eResRel = (eDep > m_threshold) - ? m_normDist() * std::sqrt( - std::pow(eRes[0] / std::sqrt(eDep), 2) + - std::pow(eRes[1], 2) + - std::pow(eRes[2] / (eDep), 2) - ) - : 0; - - const double ped = m_pedMeanADC + m_normDist() * m_pedSigmaADC; - const long long adc = std::llround(ped + eDep * (m_corrMeanScale + eResRel) / dyRangeADC * m_capADC); - - double time = std::numeric_limits::max(); - for (const auto& c : ahit.getContributions()) { - if (c.getTime() <= time) { - time = c.getTime(); - } - } - const long long tdc = std::llround((time + m_normDist() * tRes) * stepTDC); - - edm4eic::RawCalorimeterHit rawhit( - ahit.getCellID(), - (adc > m_capADC.value() ? m_capADC.value() : adc), - tdc - ); - rawhits->push_back(rawhit); - } - } - - void signal_sum_digi() { - const auto* const simhits = m_inputHitCollection.get(); - auto* rawhits = m_outputHitCollection.createAndPut(); - - // find the hits that belong to the same group (for merging) - std::unordered_map> merge_map; - for (const auto &ahit : *simhits) { - int64_t hid = (ahit.getCellID() & id_mask) | ref_mask; - auto it = merge_map.find(hid); - - if (it == merge_map.end()) { - merge_map[hid] = {ahit}; - } else { - it->second.push_back(ahit); - } - } - - // signal sum - for (auto &[id, hits] : merge_map) { - double edep = hits[0].getEnergy(); - double time = hits[0].getContributions(0).getTime(); - double max_edep = hits[0].getEnergy(); - // sum energy, take time from the most energetic hit - for (size_t i = 1; i < hits.size(); ++i) { - edep += hits[i].getEnergy(); - if (hits[i].getEnergy() > max_edep) { - max_edep = hits[i].getEnergy(); - for (const auto& c : hits[i].getContributions()) { - if (c.getTime() <= time) { - time = c.getTime(); - } - } - } - } - - // safety check - const double eResRel = (edep > m_threshold) - ? m_normDist() * eRes[0] / std::sqrt(edep) + - m_normDist() * eRes[1] + - m_normDist() * eRes[2] / edep - : 0; - - double ped = m_pedMeanADC + m_normDist() * m_pedSigmaADC; - unsigned long long adc = std::llround(ped + edep * (1. + eResRel) / dyRangeADC * m_capADC); - unsigned long long tdc = std::llround((time + m_normDist() * tRes) * stepTDC); - - edm4eic::RawCalorimeterHit rawhit( - id, - (adc > m_capADC.value() ? m_capADC.value() : adc), - tdc - ); - rawhits->push_back(rawhit); - } - } - }; - // NOLINTNEXTLINE(cppcoreguidelines-avoid-non-const-global-variables) - DECLARE_COMPONENT(CalorimeterHitDigi) +namespace { + using AlgoBase = Jug::Algo::Algorithm; +} + +class CalorimeterHitDigi : public AlgoBase { + +public: + CalorimeterHitDigi(const std::string& name, ISvcLocator* svcLoc) : AlgoBase(name, svcLoc) {} + + virtual StatusCode configure() { + setAlgoProp("energyResolutions", u_eRes.value()); + setAlgoProp("timeResolution", m_tRes.value()); + setAlgoProp("threshold", m_threshold.value()); + setAlgoProp("capacityADC", m_capADC.value()); + setAlgoProp("dynamicRangeADC", m_dyRangeADC.value()); + setAlgoProp("pedMeanADC", m_pedMeanADC.value()); + setAlgoProp("pedSigmaADC", m_pedSigmaADC.value()); + setAlgoProp("resolutionTDC", m_resolutionTDC.value()); + setAlgoProp("scaleResponse", m_corrMeanScale.value()); + setAlgoProp("fieldRefNumbers", u_refs.value()); + setAlgoProp("geoServiceName", m_geoSvcName.value()); + setAlgoProp("readoutClass", m_readout.value()); + return StatusCode::SUCCESS; + } + +private: + // additional smearing resolutions + Gaudi::Property> u_eRes{this, "energyResolutions", {}}; // a/sqrt(E/GeV) + b + c/(E/GeV) + Gaudi::Property m_tRes{this, "timeResolution", 0.0 * dd4hep::ns}; + // single hit energy deposition threshold + Gaudi::Property m_threshold{this, "threshold", 1. * dd4hep::keV}; + + // digitization settings + Gaudi::Property m_capADC{this, "capacityADC", 8096}; + Gaudi::Property m_dyRangeADC{this, "dynamicRangeADC", 100 * dd4hep::MeV}; + Gaudi::Property m_pedMeanADC{this, "pedestalMean", 400}; + Gaudi::Property m_pedSigmaADC{this, "pedestalSigma", 3.2}; + Gaudi::Property m_resolutionTDC{this, "resolutionTDC", 0.010 * dd4hep::ns}; + Gaudi::Property m_corrMeanScale{this, "scaleResponse", 1.0}; + + // signal sums + // @TODO: implement signal sums with timing + // field names to generate id mask, the hits will be grouped by masking the field + Gaudi::Property> u_fields{this, "signalSumFields", {}}; + // ref field ids are used for the merged hits, 0 is used if nothing provided + Gaudi::Property> u_refs{this, "fieldRefNumbers", {}}; + Gaudi::Property m_geoSvcName{this, "geoServiceName", "GeoSvc"}; + Gaudi::Property m_readout{this, "readoutClass", ""}; + +}; + +// NOLINTNEXTLINE(cppcoreguidelines-avoid-non-const-global-variables) +DECLARE_COMPONENT(CalorimeterHitDigi) } // namespace Jug::Digi diff --git a/JugReco/CMakeLists.txt b/JugReco/CMakeLists.txt index e4305607e31c41f27f32dc34fe92acc096cb87dc..748d0b2ee7652e1f19675a59c36c278b50ba34ec 100644 --- a/JugReco/CMakeLists.txt +++ b/JugReco/CMakeLists.txt @@ -11,7 +11,8 @@ gaudi_add_module(JugRecoPlugins ${JugRecoPlugins_sources} LINK Gaudi::GaudiAlgLib Gaudi::GaudiKernel - JugBase + JugBase JugAlgo + algocore algocalorimetry ROOT::Core ROOT::RIO ROOT::Tree EDM4HEP::edm4hep EDM4EIC::edm4eic diff --git a/JugReco/src/components/ClusterRecoCoG.cpp b/JugReco/src/components/ClusterRecoCoG.cpp index 88dd369c557264046847d5b7e3c18b0dc6b082ad..2839ca6789fb147a6d09624ed3b0af4fca4533f5 100644 --- a/JugReco/src/components/ClusterRecoCoG.cpp +++ b/JugReco/src/components/ClusterRecoCoG.cpp @@ -5,352 +5,44 @@ * Reconstruct the cluster with Center of Gravity method * Logarithmic weighting is used for mimicing energy deposit in transverse direction * - * Author: Chao Peng (ANL), 09/27/2020 + * Author: Sylvester Joosten, Chao Peng (ANL), 09/20/2022 */ -#include -#include -#include -#include -#include "boost/algorithm/string/join.hpp" -#include "fmt/format.h" -#include +#include +#include #include "Gaudi/Property.h" -#include "GaudiAlg/GaudiAlgorithm.h" -#include "GaudiAlg/GaudiTool.h" -#include "GaudiAlg/Transformer.h" -#include "GaudiKernel/PhysicalConstants.h" -#include "GaudiKernel/RndmGenerators.h" -#include "GaudiKernel/ToolHandle.h" - -#include "DDRec/CellIDPositionConverter.h" -#include "DDRec/Surface.h" -#include "DDRec/SurfaceManager.h" - -#include "JugBase/DataHandle.h" -#include "JugBase/IGeoSvc.h" - -// Event Model related classes -#include "edm4hep/MCParticle.h" -#include "edm4hep/SimCalorimeterHitCollection.h" -#include "edm4eic/ClusterCollection.h" -#include "edm4eic/MCRecoClusterParticleAssociationCollection.h" -#include "edm4eic/ProtoClusterCollection.h" -#include "edm4eic/vector_utils.h" - -using namespace Gaudi::Units; namespace Jug::Reco { -// weighting functions (with place holders for hit energy, total energy, one parameter and module -// type enum -static double constWeight(double /*E*/, double /*tE*/, double /*p*/, int /*type*/) { return 1.0; } -static double linearWeight(double E, double /*tE*/, double /*p*/, int /*type*/) { return E; } -static double logWeight(double E, double tE, double base, int /*type*/) { - return std::max(0., base + std::log(E / tE)); +namespace { + using AlgoBase = Jug::Algo::Algorithm; } -static const std::map> weightMethods{ - {"none", constWeight}, - {"linear", linearWeight}, - {"log", logWeight}, -}; +class ClusterRecoCoG : public AlgoBase { + +public: + ClusterRecoCoG(const std::string& name, ISvcLocator* svcLoc) : AlgoBase(name, svcLoc) {} + + virtual StatusCode configure() { + setAlgoProp("samplingFraction", m_sampFrac.value()); + setAlgoProp("logWeightBase", m_logWeightBase.value()); + setAlgoProp("energyWeight", m_energyWeight.value()); + setAlgoProp("moduleDimZName", m_moduleDimZName.value()); + setAlgoProp("enableEtaBounds", m_enableEtaBounds.value()); + return StatusCode::SUCCESS; + } -/** Clustering with center of gravity method. - * - * Reconstruct the cluster with Center of Gravity method - * Logarithmic weighting is used for mimicking energy deposit in transverse direction - * - * \ingroup reco - */ -class ClusterRecoCoG : public GaudiAlgorithm { private: Gaudi::Property m_sampFrac{this, "samplingFraction", 1.0}; Gaudi::Property m_logWeightBase{this, "logWeightBase", 3.6}; - Gaudi::Property m_depthCorrection{this, "depthCorrection", 0.0}; Gaudi::Property m_energyWeight{this, "energyWeight", "log"}; Gaudi::Property m_moduleDimZName{this, "moduleDimZName", ""}; - // Constrain the cluster position eta to be within - // the eta of the contributing hits. This is useful to avoid edge effects - // for endcaps. Gaudi::Property m_enableEtaBounds{this, "enableEtaBounds", false}; - - DataHandle m_inputProto{"inputProtoClusterCollection", Gaudi::DataHandle::Reader, this}; - DataHandle m_outputClusters{"outputClusterCollection", Gaudi::DataHandle::Writer, this}; - - // Collection for MC hits when running on MC - Gaudi::Property m_mcHits{this, "mcHits", ""}; - // Optional handle to MC hits - std::unique_ptr> m_mcHits_ptr; - - // Collection for associations when running on MC - Gaudi::Property m_outputAssociations{this, "outputAssociations", ""}; - // Optional handle to MC hits - std::unique_ptr> m_outputAssociations_ptr; - - // Pointer to the geometry service - SmartIF m_geoSvc; - double m_depthCorr{0}; - std::function weightFunc; - -public: - ClusterRecoCoG(const std::string& name, ISvcLocator* svcLoc) : GaudiAlgorithm(name, svcLoc) { - declareProperty("inputProtoClusterCollection", m_inputProto, ""); - declareProperty("outputClusterCollection", m_outputClusters, ""); - } - - StatusCode initialize() override { - if (GaudiAlgorithm::initialize().isFailure()) { - return StatusCode::FAILURE; - } - - // Initialize the optional MC input hit collection if requested - if (m_mcHits != "") { - m_mcHits_ptr = - std::make_unique>(m_mcHits, Gaudi::DataHandle::Reader, - this); - } - - // Initialize the optional association collection if requested - if (m_outputAssociations != "") { - m_outputAssociations_ptr = - std::make_unique>(m_outputAssociations, Gaudi::DataHandle::Writer, - this); - } - - m_geoSvc = service("GeoSvc"); - if (!m_geoSvc) { - error() << "Unable to locate Geometry Service. " - << "Make sure you have GeoSvc and SimSvc in the right order in the configuration." << endmsg; - return StatusCode::FAILURE; - } - // update depth correction if a name is provided - if (!m_moduleDimZName.value().empty()) { - m_depthCorrection = m_geoSvc->detector()->constantAsDouble(m_moduleDimZName); - } - - // select weighting method - std::string ew = m_energyWeight.value(); - // make it case-insensitive - std::transform(ew.begin(), ew.end(), ew.begin(), [](char s) { return std::tolower(s); }); - auto it = weightMethods.find(ew); - if (it == weightMethods.end()) { - error() << fmt::format("Cannot find energy weighting method {}, choose one from [{}]", m_energyWeight, - boost::algorithm::join(weightMethods | boost::adaptors::map_keys, ", ")) - << endmsg; - return StatusCode::FAILURE; - } - weightFunc = it->second; - // info() << "z_length " << depth << endmsg; - return StatusCode::SUCCESS; - } - - StatusCode execute() override { - // input collections - const auto& proto = *m_inputProto.get(); - auto& clusters = *m_outputClusters.createAndPut(); - - // Optional input MC data - const edm4hep::SimCalorimeterHitCollection* mchits = nullptr; - if (m_mcHits_ptr) { - mchits = m_mcHits_ptr->get(); - } - - // Optional output associations - edm4eic::MCRecoClusterParticleAssociationCollection* associations = nullptr; - if (m_outputAssociations_ptr) { - associations = m_outputAssociations_ptr->createAndPut(); - } - - for (const auto& pcl : proto) { - auto cl = reconstruct(pcl); - - if (msgLevel(MSG::DEBUG)) { - debug() << cl.getNhits() << " hits: " << cl.getEnergy() / GeV << " GeV, (" << cl.getPosition().x / mm << ", " - << cl.getPosition().y / mm << ", " << cl.getPosition().z / mm << ")" << endmsg; - } - clusters.push_back(cl); - - // If mcHits are available, associate cluster with MCParticle - // 1. find proto-cluster hit with largest energy deposition - // 2. find first mchit with same CellID - // 3. assign mchit's MCParticle as cluster truth - if (m_mcHits_ptr.get() != nullptr && m_outputAssociations_ptr.get() != nullptr) { - - // 1. find pclhit with largest energy deposition - auto pclhits = pcl.getHits(); - auto pclhit = std::max_element( - pclhits.begin(), - pclhits.end(), - [](const auto& pclhit1, const auto& pclhit2) { - return pclhit1.getEnergy() < pclhit2.getEnergy(); - } - ); - - // 2. find mchit with same CellID - // find_if not working, https://github.com/AIDASoft/podio/pull/273 - //auto mchit = std::find_if( - // mchits.begin(), - // mchits.end(), - // [&pclhit](const auto& mchit1) { - // return mchit1.getCellID() == pclhit->getCellID(); - // } - //); - auto mchit = mchits->begin(); - for ( ; mchit != mchits->end(); ++mchit) { - // break loop when CellID match found - if (mchit->getCellID() == pclhit->getCellID()) { - break; - } - } - if (!(mchit != mchits->end())) { - // break if no matching hit found for this CellID - warning() << "Proto-cluster has highest energy in CellID " << pclhit->getCellID() - << ", but no mc hit with that CellID was found." << endmsg; - info() << "Proto-cluster hits: " << endmsg; - for (const auto& pclhit1: pclhits) { - info() << pclhit1.getCellID() << ": " << pclhit1.getEnergy() << endmsg; - } - info() << "MC hits: " << endmsg; - for (const auto& mchit1: *mchits) { - info() << mchit1.getCellID() << ": " << mchit1.getEnergy() << endmsg; - } - break; - } - - // 3. find mchit's MCParticle - const auto& mcp = mchit->getContributions(0).getParticle(); - - // debug output - if (msgLevel(MSG::DEBUG)) { - debug() << "cluster has largest energy in cellID: " << pclhit->getCellID() << endmsg; - debug() << "pcl hit with highest energy " << pclhit->getEnergy() << " at index " << pclhit->getObjectID().index << endmsg; - debug() << "corresponding mc hit energy " << mchit->getEnergy() << " at index " << mchit->getObjectID().index << endmsg; - debug() << "from MCParticle index " << mcp.getObjectID().index << ", PDG " << mcp.getPDG() << ", " << edm4eic::magnitude(mcp.getMomentum()) << endmsg; - } - - // set association - edm4eic::MutableMCRecoClusterParticleAssociation clusterassoc; - clusterassoc.setRecID(cl.getObjectID().index); - clusterassoc.setSimID(mcp.getObjectID().index); - clusterassoc.setWeight(1.0); - clusterassoc.setRec(cl); - //clusterassoc.setSim(mcp); - associations->push_back(clusterassoc); - } else { - if (msgLevel(MSG::DEBUG)) { - debug() << "No mcHitCollection was provided, so no truth association will be performed." << endmsg; - } - } - } - - return StatusCode::SUCCESS; - } - -private: - edm4eic::MutableCluster reconstruct(const edm4eic::ProtoCluster& pcl) const { - edm4eic::MutableCluster cl; - cl.setNhits(pcl.hits_size()); - - // no hits - if (msgLevel(MSG::DEBUG)) { - debug() << "hit size = " << pcl.hits_size() << endmsg; - } - if (pcl.hits_size() == 0) { - return cl; - } - - // calculate total energy, find the cell with the maximum energy deposit - float totalE = 0.; - float maxE = 0.; - // Used to optionally constrain the cluster eta to those of the contributing hits - float minHitEta = std::numeric_limits::max(); - float maxHitEta = std::numeric_limits::min(); - auto time = pcl.getHits()[0].getTime(); - auto timeError = pcl.getHits()[0].getTimeError(); - for (unsigned i = 0; i < pcl.getHits().size(); ++i) { - const auto& hit = pcl.getHits()[i]; - const auto weight = pcl.getWeights()[i]; - if (msgLevel(MSG::DEBUG)) { - debug() << "hit energy = " << hit.getEnergy() << " hit weight: " << weight << endmsg; - } - auto energy = hit.getEnergy() * weight; - totalE += energy; - if (energy > maxE) { - } - const float eta = edm4eic::eta(hit.getPosition()); - if (eta < minHitEta) { - minHitEta = eta; - } - if (eta > maxHitEta) { - maxHitEta = eta; - } - } - cl.setEnergy(totalE / m_sampFrac); - cl.setEnergyError(0.); - cl.setTime(time); - cl.setTimeError(timeError); - - // center of gravity with logarithmic weighting - float tw = 0.; - auto v = cl.getPosition(); - for (unsigned i = 0; i < pcl.getHits().size(); ++i) { - const auto& hit = pcl.getHits()[i]; - const auto weight = pcl.getWeights()[i]; - float w = weightFunc(hit.getEnergy() * weight, totalE, m_logWeightBase.value(), 0); - tw += w; - v = v + (hit.getPosition() * w); - } - if (tw == 0.) { - warning() << "zero total weights encountered, you may want to adjust your weighting parameter." << endmsg; - } - cl.setPosition(v / tw); - cl.setPositionError({}); // @TODO: Covariance matrix - - // Optionally constrain the cluster to the hit eta values - if (m_enableEtaBounds) { - const bool overflow = (edm4eic::eta(cl.getPosition()) > maxHitEta); - const bool underflow = (edm4eic::eta(cl.getPosition()) < minHitEta); - if (overflow || underflow) { - const double newEta = overflow ? maxHitEta : minHitEta; - const double newTheta = edm4eic::etaToAngle(newEta); - const double newR = edm4eic::magnitude(cl.getPosition()); - const double newPhi = edm4eic::angleAzimuthal(cl.getPosition()); - cl.setPosition(edm4eic::sphericalToVector(newR, newTheta, newPhi)); - if (msgLevel(MSG::DEBUG)) { - debug() << "Bound cluster position to contributing hits due to " << (overflow ? "overflow" : "underflow") - << endmsg; - } - } - } - - // Additional convenience variables - - // best estimate on the cluster direction is the cluster position - // for simple 2D CoG clustering - cl.setIntrinsicTheta(edm4eic::anglePolar(cl.getPosition())); - cl.setIntrinsicPhi(edm4eic::angleAzimuthal(cl.getPosition())); - // TODO errors - - // Calculate radius - // @TODO: add skewness - if (cl.getNhits() > 1) { - double radius = 0; - for (const auto& hit : pcl.getHits()) { - const auto delta = cl.getPosition() - hit.getPosition(); - radius += delta * delta; - } - radius = sqrt((1. / (cl.getNhits() - 1.)) * radius); - cl.addToShapeParameters(radius); - cl.addToShapeParameters(0 /* skewness */); // skewness not yet calculated - } - - return cl; - } }; // NOLINTNEXTLINE(cppcoreguidelines-avoid-non-const-global-variables) DECLARE_COMPONENT(ClusterRecoCoG) } // namespace Jug::Reco + diff --git a/external/algorithms/.clang-format b/external/algorithms/.clang-format new file mode 100644 index 0000000000000000000000000000000000000000..5c443be69c45337924324a1f56fad9721dac62f6 --- /dev/null +++ b/external/algorithms/.clang-format @@ -0,0 +1,15 @@ +--- +BasedOnStyle: LLVM +BreakConstructorInitializersBeforeComma: true +ConstructorInitializerAllOnOneLineOrOnePerLine: true +Cpp11BracedListStyle: true +Standard: Cpp11 +#SpaceBeforeParens: ControlStatements +SpaceAfterControlStatementKeyword: true +PointerBindsToType: true +IncludeBlocks: Preserve +UseTab: Never +ColumnLimit: 100 +NamespaceIndentation: Inner +AlignConsecutiveAssignments: true +... diff --git a/external/algorithms/CMakeLists.txt b/external/algorithms/CMakeLists.txt new file mode 100644 index 0000000000000000000000000000000000000000..bb85179db36b93722af53cd6775e99f58ac46b40 --- /dev/null +++ b/external/algorithms/CMakeLists.txt @@ -0,0 +1,63 @@ +# SPDX-License-Identifier: LGPL-3.0-or-later +# Copyright (C) 2022 Whitney Armstrong, Wouter Deconinck, Sylvester Joosten + +cmake_minimum_required(VERSION 3.19) + +# CMP0074: find_package() uses _ROOT variables +cmake_policy(SET CMP0074 NEW) + +project(algorithms VERSION 1.0.0) + +set(CMAKE_CXX_STANDARD 17 CACHE STRING "") +if(NOT CMAKE_CXX_STANDARD MATCHES "17|20") + message(FATAL_ERROR "Unsupported C++ standard: ${CMAKE_CXX_STANDARD}") +endif() +set(CMAKE_CXX_STANDARD_REQUIRED ON) +set(CMAKE_CXX_EXTENSIONS OFF) + +# Export compile commands as json for run-clang-tidy +set(CMAKE_EXPORT_COMPILE_COMMANDS ON) + +# Also use clang-tidy integration in CMake +option(ENABLE_CLANG_TIDY "Enable clang-tidy integration in cmake" OFF) +if(ENABLE_CLANG_TIDY) + find_program(CLANG_TIDY_EXE NAMES "clang-tidy") + if (CLANG_TIDY_EXE) + message(STATUS "clang-tidy found: ${CLANG_TIDY_EXE}") + set(CMAKE_CXX_CLANG_TIDY "${CLANG_TIDY_EXE}" CACHE STRING "" FORCE) + else() + set(CMAKE_CXX_CLANG_TIDY "" CACHE STRING "" FORCE) + endif() +endif() + +# Set default build type +set(default_build_type "Release") +if(EXISTS "${CMAKE_SOURCE_DIR}/.git") + set(default_build_type "RelWithDebInfo") +endif() +if(NOT CMAKE_BUILD_TYPE AND NOT CMAKE_CONFIGURATION_TYPES) + message(STATUS "Setting build type to '${default_build_type}' as none was specified.") + set(CMAKE_BUILD_TYPE "${default_build_type}" CACHE + STRING "Choose the type of build." FORCE) + # Set the possible values of build type for cmake-gui + set_property(CACHE CMAKE_BUILD_TYPE PROPERTY STRINGS + "Debug" "Release" "MinSizeRel" "RelWithDebInfo") +endif() + +# Set all warnings +if(NOT CMAKE_BUILD_TYPE MATCHES Release) + add_compile_options(-Wall -Wextra -Werror) +endif() + +find_package(Microsoft.GSL CONFIG) + +find_package(EDM4EIC REQUIRED) +find_package(EDM4HEP 0.4.1 REQUIRED) +#find_package(DD4hep COMPONENTS DDG4 DDG4IO DDRec REQUIRED) +find_package(DD4hep COMPONENTS DDRec REQUIRED) +find_package(fmt REQUIRED) + +include(GNUInstallDirs) + +add_subdirectory(core) +add_subdirectory(calorimetry) diff --git a/external/algorithms/calorimetry/CMakeLists.txt b/external/algorithms/calorimetry/CMakeLists.txt new file mode 100644 index 0000000000000000000000000000000000000000..63890cc286a70f1f46ccdedd67be7aaa6f68d298 --- /dev/null +++ b/external/algorithms/calorimetry/CMakeLists.txt @@ -0,0 +1,45 @@ +# SPDX-License-Identifier: LGPL-3.0-or-later +# Copyright (C) 2022 Sylvester Joosten + +################################################################################ +# Package: algorithms core utilities +################################################################################ + +set(SUBDIR "calorimetry") +set(LIBRARY "algo${SUBDIR}") +set(TARGETS ${TARGETS} ${LIBRARY} PARENT_SCOPE) + +file(GLOB SRC CONFIGURE_DEPENDS src/*.cpp) + +add_library(${LIBRARY} SHARED ${SRC}) +target_link_libraries(${LIBRARY} + PUBLIC + EDM4HEP::edm4hep + EDM4EIC::edm4eic + DD4hep::DDRec + algocore + fmt::fmt) +target_include_directories(${LIBRARY} + PUBLIC + $ + $) +set_target_properties(${LIBRARY} PROPERTIES + VERSION ${PROJECT_VERSION} + SOVERSION ${PROJECT_VERSION_MAJOR}.${PROJECT_VERSION_MINOR}) + +install(TARGETS ${LIBRARY} + EXPORT algorithmsTargets + RUNTIME DESTINATION "${CMAKE_INSTALL_BINDIR}" COMPONENT bin + LIBRARY DESTINATION "${CMAKE_INSTALL_LIBDIR}" COMPONENT shlib + ARCHIVE DESTINATION "${CMAKE_INSTALL_LIBDIR}" COMPONENT lib + INCLUDES DESTINATION "${CMAKE_INSTALL_INCLUDEDIR}" + NAMESPACE algorithms::) + +install(DIRECTORY ${PROJECT_SOURCE_DIR}/${SUBDIR}/include/algorithms +DESTINATION ${CMAKE_INSTALL_INCLUDEDIR} COMPONENT dev) + +# TODO: Testing +#if(BUILD_TESTING) +# enable_testing() +#endif() + diff --git a/external/algorithms/calorimetry/include/algorithms/calorimetry/CalorimeterHitDigi.h b/external/algorithms/calorimetry/include/algorithms/calorimetry/CalorimeterHitDigi.h new file mode 100644 index 0000000000000000000000000000000000000000..c86431f9f78693371bdf1df225d8f077ab391b59 --- /dev/null +++ b/external/algorithms/calorimetry/include/algorithms/calorimetry/CalorimeterHitDigi.h @@ -0,0 +1,90 @@ +// SPDX-License-Identifier: LGPL-3.0-or-later +// Copyright (C) 2022 Chao Peng, Wouter Deconinck, Sylvester Joosten + +// A general digitization for CalorimeterHit from simulation +// 1. Smear energy deposit with a/sqrt(E/GeV) + b + c/E or a/sqrt(E/GeV) (relative value) +// 2. Digitize the energy with dynamic ADC range and add pedestal (mean +- sigma) +// 3. Time conversion with smearing resolution (absolute value) +// 4. Signal is summed if the SumFields are provided +// +// Author: Chao Peng +// Date: 06/02/2021 + +#include +#include +#include +#include + +#include "DDRec/CellIDPositionConverter.h" +#include "DDSegmentation/BitFieldCoder.h" + +#include "fmt/format.h" +#include "fmt/ranges.h" + +// Event Model related classes +#include "edm4hep/SimCalorimeterHitCollection.h" +#include "edm4eic/RawCalorimeterHitCollection.h" + +namespace algorithms::calorimetry { + + using CalorimeterHitDigiAlgorithm = Algorithm< + Input, + Output>; + + /** Generic calorimeter hit digitiziation. + * + * \ingroup digi + * \ingroup calorimetry + */ + class CalorimeterHitDigi : public CalorimeterHitDigiAlgorithm { + public: + using Input = CalorimeterHitDigiAlgorithm::Input; + using Output = CalorimeterHitDigiAlgorithm::Output; + + CalorimeterHitDigi(std::string_view name) + : CalorimeterHitDigiAlgorithm{name, + {"inputHitCollection"}, + {"outputHitCollection"}} {}; + + void init(); + void process(const Input&, const Output&); + + private: + void single_hits_digi(const Input&, const Output&); + void signal_sum_digi(const Input&, const Output&); + + // additional smearing resolutions + Property> u_eRes{this, "energyResolutions", {}}; // a/sqrt(E/GeV) + b + c/(E/GeV) + Property m_tRes{this, "timeResolution", 0.0 * dd4hep::ns}; + // single hit energy deposition threshold + Property m_threshold{this, "threshold", 1. * dd4hep::keV}; + + // digitization settings + Property m_capADC{this, "capacityADC", 8096}; + Property m_dyRangeADC{this, "dynamicRangeADC", 100 * dd4hep::MeV}; + Property m_pedMeanADC{this, "pedestalMean", 400}; + Property m_pedSigmaADC{this, "pedestalSigma", 3.2}; + Property m_resolutionTDC{this, "resolutionTDC", 0.010 * dd4hep::ns}; + + Property m_corrMeanScale{this, "scaleResponse", 1.0}; + // These are better variable names for the "energyResolutions" array which is a bit + // magic @FIXME + //Property m_corrSigmaCoeffE{this, "responseCorrectionSigmaCoeffE", 0.0}; + //Property m_corrSigmaCoeffSqrtE{this, "responseCorrectionSigmaCoeffSqrtE", 0.0}; + + // signal sums + // @TODO: implement signal sums with timing + // field names to generate id mask, the hits will be grouped by masking the field + Property> u_fields{this, "signalSumFields", {}}; + // ref field ids are used for the merged hits, 0 is used if nothing provided + Property> u_refs{this, "fieldRefNumbers", {}}; + Property m_readout{this, "readoutClass", ""}; + + double eRes[3] = {0., 0., 0.}; + uint64_t id_mask{0}, ref_mask{0}; + + const GeoSvc& m_geoSvc = GeoSvc::instance(); + const RandomSvc& m_randomSvc = RandomSvc::instance(); + }; + +} // namespace algoriths::calorimetry diff --git a/external/algorithms/calorimetry/include/algorithms/calorimetry/ClusterRecoCoG.h b/external/algorithms/calorimetry/include/algorithms/calorimetry/ClusterRecoCoG.h new file mode 100644 index 0000000000000000000000000000000000000000..00afe78470aaafda0a704cf30a52b776ed14992a --- /dev/null +++ b/external/algorithms/calorimetry/include/algorithms/calorimetry/ClusterRecoCoG.h @@ -0,0 +1,67 @@ +// SPDX-License-Identifier: LGPL-3.0-or-later +// Copyright (C) 2022 Sylvester Joosten, Chao, Chao Peng, Whitney Armstrong + +/* + * Reconstruct the cluster with Center of Gravity method + * Logarithmic weighting is used for mimicing energy deposit in transverse direction + * + * Author: Sylvester Joosten (ANL), Chao Peng (ANL) 09/19/2022 + */ + +#include +#include +#include + +// Data types +#include +#include +#include +#include + +namespace algorithms::calorimetry { + +using ClusteringAlgorithm = Algorithm< + Input>, + Output>>; + +/** Clustering with center of gravity method. + * + * Reconstruct the cluster with Center of Gravity method + * Logarithmic weighting is used for mimicking energy deposit in transverse direction + * + * \ingroup reco + */ +class ClusterRecoCoG : public ClusteringAlgorithm { +public: + using Input = ClusteringAlgorithm::Input; + using Output = ClusteringAlgorithm::Output; + using WeightFunc = std::function; + + // TODO: get rid of "Collection" in names + ClusterRecoCoG(std::string_view name) + : ClusteringAlgorithm{name, + {"inputProtoClusterCollection", "mcHits"}, + {"outputClusterCollection", "outputAssociations"}} {} + + void init(); + void process(const Input&, const Output&); + +private: + edm4eic::MutableCluster reconstruct(const edm4eic::ProtoCluster&) const; + + Property m_sampFrac{this, "samplingFraction", 1.0}; + Property m_logWeightBase{this, "logWeightBase", 3.6}; + Property m_energyWeight{this, "energyWeight", "log"}; + Property m_moduleDimZName{this, "moduleDimZName", ""}; + // Constrain the cluster position eta to be within + // the eta of the contributing hits. This is useful to avoid edge effects + // for endcaps. + Property m_enableEtaBounds{this, "enableEtaBounds", true}; + + WeightFunc m_weightFunc; + + const GeoSvc& m_geo = GeoSvc::instance(); +}; +} // namespace algorithms::calorimetry + diff --git a/external/algorithms/calorimetry/src/CalorimeterHitDigi.cpp b/external/algorithms/calorimetry/src/CalorimeterHitDigi.cpp new file mode 100644 index 0000000000000000000000000000000000000000..deae3f9590260382acac888bc9c1aec3826fcaf0 --- /dev/null +++ b/external/algorithms/calorimetry/src/CalorimeterHitDigi.cpp @@ -0,0 +1,175 @@ +// SPDX-License-Identifier: LGPL-3.0-or-later +// Copyright (C) 2022 Chao Peng, Wouter Deconinck, Sylvester Joosten + +// A general digitization for CalorimeterHit from simulation +// 1. Smear energy deposit with a/sqrt(E/GeV) + b + c/E or a/sqrt(E/GeV) (relative value) +// 2. Digitize the energy with dynamic ADC range and add pedestal (mean +- sigma) +// 3. Time conversion with smearing resolution (absolute value) +// 4. Signal is summed if the SumFields are provided +// +// Author: Chao Peng +// Date: 06/02/2021 + +#include + +#include +#include +#include + +#include "fmt/format.h" +#include "fmt/ranges.h" + +// Event Model related classes +#include "edm4hep/SimCalorimeterHitCollection.h" +#include "edm4eic/RawCalorimeterHitCollection.h" +#include "edm4eic/RawCalorimeterHitData.h" + +namespace algorithms::calorimetry { + +void CalorimeterHitDigi::init() { + // set energy resolution numbers + for (size_t i = 0; i < u_eRes.size() && i < 3; ++i) { + eRes[i] = u_eRes[i]; + } + + // need signal sum + if (!u_fields.value().empty()) { + // sanity checks + if (m_readout.value().empty()) { + error() << "readoutClass is not provided, it is needed to know the fields in readout ids" + << endmsg; + return; + } + + // get decoders + try { + auto id_desc = m_geoSvc.detector()->readout(m_readout).idSpec(); + id_mask = 0; + std::vector> ref_fields; + for (size_t i = 0; i < u_fields.size(); ++i) { + id_mask |= id_desc.field(u_fields[i])->mask(); + // use the provided id number to find ref cell, or use 0 + int ref = i < u_refs.size() ? u_refs[i] : 0; + ref_fields.emplace_back(u_fields[i], ref); + } + ref_mask = id_desc.encode(ref_fields); + // debug() << fmt::format("Referece id mask for the fields {:#064b}", ref_mask) << endmsg; + } catch (...) { + error() << "Failed to load ID decoder for " << m_readout << endmsg; + return; + } + id_mask = ~id_mask; + info() << fmt::format("ID mask in {:s}: {:#064b}", m_readout, id_mask) << endmsg; + } +} + +void CalorimeterHitDigi::process( + const CalorimeterHitDigi::Input& input, + const CalorimeterHitDigi::Output& output +) { + if (!u_fields.value().empty()) { + signal_sum_digi(input, output); + } else { + single_hits_digi(input, output); + } +} + +void CalorimeterHitDigi::single_hits_digi( + const CalorimeterHitDigi::Input& input, + const CalorimeterHitDigi::Output& output +) { + const auto [simhits] = input; + auto [rawhits] = output; + + for (const auto& ahit : *simhits) { + // Note: juggler internal unit of energy is GeV + const double eDep = ahit.getEnergy(); + + // apply additional calorimeter noise to corrected energy deposit + const double eResRel = (eDep > m_threshold) + ? m_randomSvc.normal() * std::sqrt( + std::pow(eRes[0] / std::sqrt(eDep), 2) + + std::pow(eRes[1], 2) + + std::pow(eRes[2] / (eDep), 2) + ) + : 0; + + const double ped = m_pedMeanADC + m_randomSvc.normal() * m_pedSigmaADC; + const long long adc = std::llround(ped + eDep * (m_corrMeanScale + eResRel) / m_dyRangeADC * m_capADC); + + double time = std::numeric_limits::max(); + for (const auto& c : ahit.getContributions()) { + if (c.getTime() <= time) { + time = c.getTime(); + } + } + const long long tdc = std::llround((time + m_randomSvc.normal() * m_tRes) / m_resolutionTDC); + + edm4eic::RawCalorimeterHit rawhit( + ahit.getCellID(), + (adc > m_capADC.value() ? m_capADC.value() : adc), + tdc + ); + rawhits->push_back(rawhit); + } +} + +void CalorimeterHitDigi::signal_sum_digi( + const CalorimeterHitDigi::Input& input, + const CalorimeterHitDigi::Output& output +) { + const auto [simhits] = input; + auto [rawhits] = output; + + // find the hits that belong to the same group (for merging) + std::unordered_map> merge_map; + for (const auto &ahit : *simhits) { + int64_t hid = (ahit.getCellID() & id_mask) | ref_mask; + auto it = merge_map.find(hid); + + if (it == merge_map.end()) { + merge_map[hid] = {ahit}; + } else { + it->second.push_back(ahit); + } + } + + // signal sum + for (auto &[id, hits] : merge_map) { + double edep = hits[0].getEnergy(); + double time = hits[0].getContributions(0).getTime(); + double max_edep = hits[0].getEnergy(); + // sum energy, take time from the most energetic hit + for (size_t i = 1; i < hits.size(); ++i) { + edep += hits[i].getEnergy(); + if (hits[i].getEnergy() > max_edep) { + max_edep = hits[i].getEnergy(); + for (const auto& c : hits[i].getContributions()) { + if (c.getTime() <= time) { + time = c.getTime(); + } + } + } + } + + // safety check + const double eResRel = (edep > m_threshold) + ? m_randomSvc.normal() * eRes[0] / std::sqrt(edep) + + m_randomSvc.normal() * eRes[1] + + m_randomSvc.normal() * eRes[2] / edep + : 0; + + double ped = m_pedMeanADC + m_randomSvc.normal() * m_pedSigmaADC; + unsigned long long adc = std::llround(ped + edep * (1. + eResRel) / m_dyRangeADC * m_capADC); + unsigned long long tdc = std::llround((time + m_randomSvc.normal() * m_tRes) / m_resolutionTDC); + + edm4eic::RawCalorimeterHit rawhit( + id, + (adc > m_capADC.value() ? m_capADC.value() : adc), + tdc + ); + rawhits->push_back(rawhit); + } +} + +} // namespace algorithms::calorimetry diff --git a/external/algorithms/calorimetry/src/ClusterRecoCoG.cpp b/external/algorithms/calorimetry/src/ClusterRecoCoG.cpp new file mode 100644 index 0000000000000000000000000000000000000000..3feb2476a9841fb7a610142f248ada530d932504 --- /dev/null +++ b/external/algorithms/calorimetry/src/ClusterRecoCoG.cpp @@ -0,0 +1,248 @@ +// SPDX-License-Identifier: LGPL-3.0-or-later +// Copyright (C) 2022 Sylvester Joosten, Chao, Chao Peng, Whitney Armstrong + +/* + * Reconstruct the cluster with Center of Gravity method + * Logarithmic weighting is used for mimicing energy deposit in transverse direction + * + * Author: Chao Peng (ANL), 09/27/2020 + */ + +#include + +#include +#include +#include + +#include +#include + +// Event Model related classes +#include "edm4eic/vector_utils.h" + +namespace algorithms::calorimetry { +namespace { + + // weighting functions (with place holders for hit energy, total energy, one parameter + double constWeight(double /*E*/, double /*tE*/, double /*p*/) { return 1.0; } + double linearWeight(double E, double /*tE*/, double /*p*/) { return E; } + double logWeight(double E, double tE, double base) { + return std::max(0., base + std::log(E / tE)); + } + + const std::map weightMethods{ + {"none", constWeight}, + {"linear", linearWeight}, + {"log", logWeight}, + }; +} // namespace + +void ClusterRecoCoG::init() { + // select weighting method + std::string ew = m_energyWeight; + // make it case-insensitive + std::transform(ew.begin(), ew.end(), ew.begin(), [](char s) { return std::tolower(s); }); + if (weightMethods.count(ew)) { + std::vector keys; + std::transform(weightMethods.begin(), weightMethods.end(), std::back_inserter(keys), + [](const auto& keyvalue) { return keyvalue.first; }); + raise(fmt::format("Cannot find energy weighting method {}, choose one from {}", m_energyWeight, + keys)); + } + m_weightFunc = weightMethods.at(ew); + info() << fmt::format("Energy weight method set to: {}", ew) << endmsg; +} + +void ClusterRecoCoG::process(const ClusterRecoCoG::Input& input, + const ClusterRecoCoG::Output& output) { + const auto [proto, opt_simhits] = input; + auto [clusters, opt_assoc] = output; + + for (const auto& pcl : *proto) { + auto cl = reconstruct(pcl); + + if (aboveDebugThreshold()) { + debug() << cl.getNhits() << " hits: " << cl.getEnergy() / dd4hep::GeV << " GeV, (" + << cl.getPosition().x / dd4hep::mm << ", " << cl.getPosition().y / dd4hep::mm << ", " + << cl.getPosition().z / dd4hep::mm << ")" << endmsg; + } + clusters->push_back(cl); + + // If mcHits are available, associate cluster with MCParticle + // 1. find proto-cluster hit with largest energy deposition + // 2. find first mchit with same CellID + // 3. assign mchit's MCParticle as cluster truth + if (opt_simhits && opt_assoc) { + + // 1. find pclhit with largest energy deposition + auto pclhits = pcl.getHits(); + auto pclhit = std::max_element(pclhits.begin(), pclhits.end(), + [](const auto& pclhit1, const auto& pclhit2) { + return pclhit1.getEnergy() < pclhit2.getEnergy(); + }); + + // 2. find mchit with same CellID + // find_if not working, https://github.com/AIDASoft/podio/pull/273 + // auto mchit = std::find_if( + // opt_simhits->begin(), + // opt_simhits->end(), + // [&pclhit](const auto& mchit1) { + // return mchit1.getCellID() == pclhit->getCellID(); + // } + //); + auto mchit = opt_simhits->begin(); + for (; mchit != opt_simhits->end(); ++mchit) { + // break loop when CellID match found + if (mchit->getCellID() == pclhit->getCellID()) { + break; + } + } + if (!(mchit != opt_simhits->end())) { + // error condition should not happen + // break if no matching hit found for this CellID + warning() << "Proto-cluster has highest energy in CellID " << pclhit->getCellID() + << ", but no mc hit with that CellID was found." << endmsg; + info() << "Proto-cluster hits: " << endmsg; + for (const auto& pclhit1 : pclhits) { + info() << pclhit1.getCellID() << ": " << pclhit1.getEnergy() << endmsg; + } + info() << "MC hits: " << endmsg; + for (const auto& mchit1 : *opt_simhits) { + info() << mchit1.getCellID() << ": " << mchit1.getEnergy() << endmsg; + } + break; + } + + // 3. find mchit's MCParticle + const auto& mcp = mchit->getContributions(0).getParticle(); + + // debug output + if (aboveDebugThreshold()) { + debug() << "cluster has largest energy in cellID: " << pclhit->getCellID() << endmsg; + debug() << "pcl hit with highest energy " << pclhit->getEnergy() << " at index " + << pclhit->getObjectID().index << endmsg; + debug() << "corresponding mc hit energy " << mchit->getEnergy() << " at index " + << mchit->getObjectID().index << endmsg; + debug() << "from MCParticle index " << mcp.getObjectID().index << ", PDG " << mcp.getPDG() + << ", " << edm4eic::magnitude(mcp.getMomentum()) << endmsg; + } + + // set association + edm4eic::MutableMCRecoClusterParticleAssociation clusterassoc; + clusterassoc.setRecID(cl.getObjectID().index); + clusterassoc.setSimID(mcp.getObjectID().index); + clusterassoc.setWeight(1.0); + clusterassoc.setRec(cl); + // clusterassoc.setSim(mcp); + opt_assoc->push_back(clusterassoc); + } else { + if (aboveDebugThreshold()) { + debug() << "No mcHitCollection was provided, so no truth association will be performed." + << endmsg; + } + } + } +} + +edm4eic::MutableCluster ClusterRecoCoG::reconstruct(const edm4eic::ProtoCluster& pcl) const { + edm4eic::MutableCluster cl; + cl.setNhits(pcl.hits_size()); + + // no hits + if (aboveDebugThreshold()) { + debug() << "hit size = " << pcl.hits_size() << endmsg; + } + if (pcl.hits_size() == 0) { + return cl; + } + + // calculate total energy, find the cell with the maximum energy deposit + float totalE = 0.; + float maxE = 0.; + // Used to optionally constrain the cluster eta to those of the contributing hits + float minHitEta = std::numeric_limits::max(); + float maxHitEta = std::numeric_limits::min(); + auto time = pcl.getHits()[0].getTime(); + auto timeError = pcl.getHits()[0].getTimeError(); + for (unsigned i = 0; i < pcl.getHits().size(); ++i) { + const auto& hit = pcl.getHits()[i]; + const auto weight = pcl.getWeights()[i]; + if (aboveDebugThreshold()) { + debug() << "hit energy = " << hit.getEnergy() << " hit weight: " << weight << endmsg; + } + auto energy = hit.getEnergy() * weight; + totalE += energy; + if (energy > maxE) { + } + const float eta = edm4eic::eta(hit.getPosition()); + if (eta < minHitEta) { + minHitEta = eta; + } + if (eta > maxHitEta) { + maxHitEta = eta; + } + } + cl.setEnergy(totalE / m_sampFrac); + cl.setEnergyError(0.); + cl.setTime(time); + cl.setTimeError(timeError); + + // center of gravity with logarithmic weighting + float tw = 0.; + auto v = cl.getPosition(); + for (unsigned i = 0; i < pcl.getHits().size(); ++i) { + const auto& hit = pcl.getHits()[i]; + const auto weight = pcl.getWeights()[i]; + float w = m_weightFunc(hit.getEnergy() * weight, totalE, m_logWeightBase.value()); + tw += w; + v = v + (hit.getPosition() * w); + } + if (tw == 0.) { + warning() << "zero total weights encountered, you may want to adjust your weighting parameter." + << endmsg; + } + cl.setPosition(v / tw); + cl.setPositionError({}); // @TODO: Covariance matrix + + // Optionally constrain the cluster to the hit eta values + if (m_enableEtaBounds) { + const bool overflow = (edm4eic::eta(cl.getPosition()) > maxHitEta); + const bool underflow = (edm4eic::eta(cl.getPosition()) < minHitEta); + if (overflow || underflow) { + const double newEta = overflow ? maxHitEta : minHitEta; + const double newTheta = edm4eic::etaToAngle(newEta); + const double newR = edm4eic::magnitude(cl.getPosition()); + const double newPhi = edm4eic::angleAzimuthal(cl.getPosition()); + cl.setPosition(edm4eic::sphericalToVector(newR, newTheta, newPhi)); + if (aboveDebugThreshold()) { + debug() << "Bound cluster position to contributing hits due to " + << (overflow ? "overflow" : "underflow") << endmsg; + } + } + } + + // Additional convenience variables + + // best estimate on the cluster direction is the cluster position + // for simple 2D CoG clustering + cl.setIntrinsicTheta(edm4eic::anglePolar(cl.getPosition())); + cl.setIntrinsicPhi(edm4eic::angleAzimuthal(cl.getPosition())); + // TODO errors + + // Calculate radius + // @TODO: add skewness + if (cl.getNhits() > 1) { + double radius = 0; + for (const auto& hit : pcl.getHits()) { + const auto delta = cl.getPosition() - hit.getPosition(); + radius += delta * delta; + } + radius = sqrt((1. / (cl.getNhits() - 1.)) * radius); + cl.addToShapeParameters(radius); + cl.addToShapeParameters(0 /* skewness */); // skewness not yet calculated + } + + return cl; +} + +} // namespace algorithms::calorimetry diff --git a/external/algorithms/cmake/algorithmsConfig.cmake b/external/algorithms/cmake/algorithmsConfig.cmake new file mode 100644 index 0000000000000000000000000000000000000000..2b5d16dab0ee6da7f467ee51681e551db294b556 --- /dev/null +++ b/external/algorithms/cmake/algorithmsConfig.cmake @@ -0,0 +1,5 @@ +include(CMakeFindDependencyMacro) + +# - Include the targets file to create the imported targets that a client can +# link to (libraries) or execute (programs) +include("${CMAKE_CURRENT_LIST_DIR}/algorithmsTargets.cmake") diff --git a/external/algorithms/core/CMakeLists.txt b/external/algorithms/core/CMakeLists.txt new file mode 100644 index 0000000000000000000000000000000000000000..153f9a282e53415c1f7f0d85fc5fcf0f112a29b8 --- /dev/null +++ b/external/algorithms/core/CMakeLists.txt @@ -0,0 +1,44 @@ +# SPDX-License-Identifier: LGPL-3.0-or-later +# Copyright (C) 2022 Sylvester Joosten + +################################################################################ +# Package: algorithms core utilities +################################################################################ + +set(SUBDIR "core") +set(LIBRARY "algo${SUBDIR}") +set(TARGETS ${TARGETS} ${LIBRARY} PARENT_SCOPE) + +file(GLOB SRC CONFIGURE_DEPENDS src/*.cpp) + +add_library(${LIBRARY} SHARED ${SRC}) +target_link_libraries(${LIBRARY} + PUBLIC + EDM4HEP::edm4hep + EDM4EIC::edm4eic + DD4hep::DDRec + fmt::fmt) +target_include_directories(${LIBRARY} + PUBLIC + $ + $) +set_target_properties(${LIBRARY} PROPERTIES + VERSION ${PROJECT_VERSION} + SOVERSION ${PROJECT_VERSION_MAJOR}.${PROJECT_VERSION_MINOR}) + +install(TARGETS ${LIBRARY} + EXPORT algorithmsTargets + RUNTIME DESTINATION "${CMAKE_INSTALL_BINDIR}" COMPONENT bin + LIBRARY DESTINATION "${CMAKE_INSTALL_LIBDIR}" COMPONENT shlib + ARCHIVE DESTINATION "${CMAKE_INSTALL_LIBDIR}" COMPONENT lib + INCLUDES DESTINATION "${CMAKE_INSTALL_INCLUDEDIR}" + NAMESPACE algorithms::) + +install(DIRECTORY ${PROJECT_SOURCE_DIR}/${SUBDIR}/include/algorithms +DESTINATION ${CMAKE_INSTALL_INCLUDEDIR} COMPONENT dev) + +# TODO: Testing +#if(BUILD_TESTING) +# enable_testing() +#endif() + diff --git a/external/algorithms/core/include/algorithms/algorithm.h b/external/algorithms/core/include/algorithms/algorithm.h new file mode 100644 index 0000000000000000000000000000000000000000..cad074fc634d9a108d3642cdcfb1984af1b1bf4e --- /dev/null +++ b/external/algorithms/core/include/algorithms/algorithm.h @@ -0,0 +1,69 @@ +#pragma once + +#include +#include +#include +#include +#include + +#include +#include +#include +#include + +namespace algorithms { + +// T should either be the desired input type, a std::vector<> of the desired input type, +// or a std::optional<> of the desired input type +template struct Input : std::tuple...> { + constexpr static const size_t kSize = sizeof...(T); + using value_type = std::tuple...>; + using data_type = std::tuple; + using key_type = std::array; +}; +template struct Output : std::tuple...> { + constexpr static const size_t kSize = sizeof...(T); + using value_type = std::tuple...>; + using data_type = std::tuple; + using key_type = std::array; +}; + +// TODO: C++20 Concepts version for better error handling +template +class Algorithm : public PropertyMixin, public LoggerMixin, public NameMixin { +public: + using input_type = InputType; + using output_type = OutputType; + using Input = typename input_type::value_type; + using Output = typename output_type::value_type; + using InputNames = typename input_type::key_type; + using OutputNames = typename output_type::key_type; + + Algorithm(std::string_view name, const InputNames& input_names, const OutputNames& output_names) + : LoggerMixin(name) + , NameMixin(name) + , m_input_names{input_names} + , m_output_names{output_names} {} + + void init(); + void process(const Input& input, const Output& output); + + const InputNames& inputNames() const { return m_input_names; } + const OutputNames& outputNames() const { return m_output_names; } + +private: + const InputNames m_input_names; + const OutputNames m_output_names; +}; + +namespace detail { + template struct is_input : std::false_type {}; + template struct is_input> : std::true_type {}; + template struct is_output : std::false_type {}; + template struct is_output> : std::true_type {}; +} // namespace detail +template constexpr bool is_input_v = detail::is_input::value; +template constexpr bool is_output_v = detail::is_output::value; + +} // namespace algorithms + diff --git a/external/algorithms/core/include/algorithms/detail/upcast.h b/external/algorithms/core/include/algorithms/detail/upcast.h new file mode 100644 index 0000000000000000000000000000000000000000..4cbe1c0dcb260b4e710ecc19e6b99a300e182ee5 --- /dev/null +++ b/external/algorithms/core/include/algorithms/detail/upcast.h @@ -0,0 +1,35 @@ +#pragma once + +#include +#include +#include + +// Safe type conversions (without narrowing) to be used for properties, +// gets arround the limitations of type-erasure without explicitly needing +// to specify template arguments in setProperty etc. + +namespace algorithms::detail { +template struct upcast_type { using type = T; }; +template +struct upcast_type && !std::is_signed_v>> { + using type = uint64_t; +}; +template +struct upcast_type && std::is_signed_v>> { + using type = int64_t; +}; +template +struct upcast_type>> { + using type = double; +}; +template +struct upcast_type>> { + using type = std::string; +}; + +template using upcast_type_t = typename upcast_type::type; + +template upcast_type_t upcast(T&& value) { return value; } + +} // namespace algorithms::detail + diff --git a/external/algorithms/core/include/algorithms/error.h b/external/algorithms/core/include/algorithms/error.h new file mode 100644 index 0000000000000000000000000000000000000000..195c03224833de10f539ae537f835400696cfa20 --- /dev/null +++ b/external/algorithms/core/include/algorithms/error.h @@ -0,0 +1,23 @@ +#pragma once + +#include +#include + +namespace algorithms { + +class Error : public std::exception { +public: + Error(std::string_view msg, std::string_view type = "algorithms::Error") + : m_msg{msg}, m_type{type} {} + + virtual const char* what() const noexcept { return m_msg.c_str(); } + virtual const char* type() const noexcept { return m_type.c_str(); } + virtual ~Error() noexcept {} + +private: + const std::string m_msg; + const std::string m_type; +}; + +} // namespace algorithms + diff --git a/external/algorithms/core/include/algorithms/geo.h b/external/algorithms/core/include/algorithms/geo.h new file mode 100644 index 0000000000000000000000000000000000000000..319ad531af557baa39a06c08ebd71648f22da2b0 --- /dev/null +++ b/external/algorithms/core/include/algorithms/geo.h @@ -0,0 +1,39 @@ +#pragma once + +#include +#include + +#include "DDRec/CellIDPositionConverter.h" +#include + +#include +#include + +namespace algorithms { + +class GeoSvc : public LoggedService { +public: + // Initialize the geometry service, to be called after with a global detector + // pointer (will not re-initialize), or after at least the detectors property was specified + // (will load XML and then fully initialize) + void init(dd4hep::Detector* = nullptr); + + // TODO check const-ness + gsl::not_null detector() const { return m_detector; } + dd4hep::DetElement world() const { return detector()->world(); } + gsl::not_null cellIDPositionConverter() const { + return m_converter.get(); + } + +private: + dd4hep::Detector* m_detector = nullptr; + std::unique_ptr m_converter; + + // Configuration variables. These only need to be specified if we are actually running + // in "standalone" mode (hence they're optional) + Property> m_xml_list{this, "detectors", {}}; + + ALGORITHMS_DEFINE_LOGGED_SERVICE(GeoSvc) +}; + +} // namespace algorithms diff --git a/external/algorithms/core/include/algorithms/logger.h b/external/algorithms/core/include/algorithms/logger.h new file mode 100644 index 0000000000000000000000000000000000000000..a46df4f9b635d67157ca2217f1064db5d76acaf3 --- /dev/null +++ b/external/algorithms/core/include/algorithms/logger.h @@ -0,0 +1,234 @@ +#pragma once + +#include +#include +#include +#include +#include +#include +#include +#include + +#include +#include +#include + +// Simple thread-safe logger with optional overrides by the calling framework +namespace algorithms { + +enum class LogLevel : unsigned { + kTrace = 0, + kDebug = 1, + kInfo = 2, + kWarning = 3, + kError = 4, + kCritical = 5, +}; +constexpr std::string_view logLevelName(LogLevel level) { + // Compiler can warn if not all of the enum is covered + switch (level) { + case LogLevel::kTrace: + return "TRACE"; + case LogLevel::kDebug: + return "DEBUG"; + case LogLevel::kInfo: + return "INFO"; + case LogLevel::kWarning: + return "WARNING"; + case LogLevel::kError: + return "ERROR"; + case LogLevel::kCritical: + return "CRITICAL"; + } + // Default return to make gcc happy, will never happen + return "UNKNOWN"; +} + +// Note: the log action is responsible for dealing with concurrent calls +// the default LogAction is a thread-safe example +class LogSvc : public Service { +public: + using LogAction = std::function; + void defaultLevel(const LogLevel l) { m_level.set(l); } + LogLevel defaultLevel() const { return m_level; } + void action(LogAction a) { m_action = a; } + void report(const LogLevel l, std::string_view caller, std::string_view msg) const { + m_action(l, caller, msg); + } + +private: + Property m_level{this, "defaultLevel", LogLevel::kInfo}; + LogAction m_action = [](const LogLevel l, std::string_view caller, std::string_view msg) { + static std::mutex m; + std::lock_guard lock(m); + fmt::print("{} [{}] {}\n", logLevelName(l), caller, msg); + }; + ALGORITHMS_DEFINE_SERVICE(LogSvc) +}; + +namespace detail { + // Output buffer that calls our global logger's report() function + class LoggerBuffer : public std::stringbuf { + public: + LoggerBuffer(const LogLevel l, std::string_view caller) + : m_mylevel{l}, m_caller{caller}, m_logger{LogSvc::instance()} {} + virtual int sync() { + // report should deal with concurrency (the minimal version does) + m_logger.report(m_mylevel, m_caller, this->str()); + this->str(""); + return 0; + } + + private: + // The output buffer knows the log level of its associated logger + // (eg. is this the debug logger?) + LogLevel m_mylevel; + const std::string m_caller; + const LogSvc& m_logger; + }; + // thread-safe output stream for the logger + class LoggerStream { + public: + LoggerStream(std::string_view caller, const LogLevel level, + const LogLevel threshold = LogSvc::instance().defaultLevel()) + : m_buffer{level, caller}, m_os{&m_buffer}, m_level{level}, m_threshold{threshold} {} + LoggerStream() = delete; + LoggerStream(const LoggerStream&) = delete; + + template LoggerStream& operator<<(Arg&& streamable) { + if (m_level >= m_threshold) { + std::lock_guard lock{m_mutex}; + m_os << std::forward(streamable); + return *this; + } + return *this; + } + // To support input manipulators such as std::endl + // Note: would be better with Concepts + using IOManipType1 = std::ostream&(std::ostream&); // this capturs std::endl; + using IOManipType2 = std::ios_base&(std::ios_base&); + LoggerStream& operator<<(IOManipType1* f) { + if (m_level >= m_threshold) { + std::lock_guard lock{m_mutex}; + f(m_os); + return *this; + } + return *this; + } + LoggerStream& operator<<(IOManipType2* f) { + if (m_level >= m_threshold) { + std::lock_guard lock{m_mutex}; + f(m_os); + return *this; + } + return *this; + } + LogLevel threshold() const { return m_threshold; } + void threshold(const LogLevel th) { m_threshold = th; } + + private: + std::mutex m_mutex; + LoggerBuffer m_buffer; + std::ostream m_os; + const LogLevel m_level; + LogLevel m_threshold; + }; +} // namespace detail + +// Mixin meant to add utility logger functions to algorithms/services/etc +class LoggerMixin { +public: + LoggerMixin(std::string_view caller, const LogLevel threshold = LogSvc::instance().defaultLevel()) + : m_caller{caller}, m_logger{LogSvc::instance()} { + level(threshold); + } + +public: + // Not done through Properties, as that would require entanglement with the + // PropertyMixin which is not appropriate here. It's the FW responsibility to set this + // on the algorithm level if desired, before or during the init() stage. + void level(const LogLevel threshold) { + m_level = threshold; + m_critical.threshold(m_level); + m_error.threshold(m_level); + m_warning.threshold(m_level); + m_info.threshold(m_level); + m_debug.threshold(m_level); + m_trace.threshold(m_level); + } + LogLevel level() const { return m_level; } + +protected: + detail::LoggerStream& critical() const { return m_critical; } + detail::LoggerStream& error() const { return m_error; } + detail::LoggerStream& warning() const { return m_warning; } + detail::LoggerStream& info() const { return m_info; } + detail::LoggerStream& debug() const { return m_debug; } + detail::LoggerStream& trace() const { return m_trace; } + + void critical(std::string_view msg) { report(msg); } + void error(std::string_view msg) { report(msg); } + void warning(std::string_view msg) { report(msg); } + void info(std::string_view msg) { report(msg); } + void debug(std::string_view msg) { report(msg); } + void trace(std::string_view msg) { report(msg); } + + bool aboveCriticalThreshold() const { return m_level >= LogLevel::kCritical; } + bool aboveErrorThreshold() const { return m_level >= LogLevel::kError; } + bool aboveWarningThreshold() const { return m_level >= LogLevel::kWarning; } + bool aboveInfoThreshold() const { return m_level >= LogLevel::kInfo; } + bool aboveDebugThreshold() const { return m_level >= LogLevel::kDebug; } + bool aboveTraceThreshold() const { return m_level >= LogLevel::kTrace; } + + // LoggerMixin also provides nice error raising + // ErrorTypes needs to derive from Error, and needs to have a constructor that takes two + // strings --> TODO add C++20 Concept version + template void raise(std::string_view msg) const { + error() << msg << endmsg; + throw ErrorType(msg); + } + + // Endmsg (flush) support to initiate a sync() action + static std::ostream& endmsg(std::ostream& os) { + os.flush(); + return os; + } + +private: + template void report(std::string_view msg) { + if (l >= m_level) { + m_logger.report(l, m_caller, msg); + } + } + + const std::string m_caller; + LogLevel m_level; + mutable detail::LoggerStream m_critical{m_caller, LogLevel::kCritical}; + mutable detail::LoggerStream m_error{m_caller, LogLevel::kError}; + mutable detail::LoggerStream m_warning{m_caller, LogLevel::kWarning}; + mutable detail::LoggerStream m_info{m_caller, LogLevel::kInfo}; + mutable detail::LoggerStream m_debug{m_caller, LogLevel::kDebug}; + mutable detail::LoggerStream m_trace{m_caller, LogLevel::kTrace}; + + const LogSvc& m_logger; +}; + +// A Service base class with logger functionality +template class LoggedService : public Service, public LoggerMixin { +protected: + LoggedService(std::string_view name) : Service(name), LoggerMixin(name) {} +}; + +} // namespace algorithms + +#define ALGORITHMS_DEFINE_LOGGED_SERVICE(className) \ +protected: \ + className() : LoggedService(#className) {} \ + \ +public: \ + friend class Service; \ + className(const className&) = delete; \ + void operator=(const className&) = delete; \ + constexpr static const char* kName = #className; + +//#define endmsg std::flush diff --git a/external/algorithms/core/include/algorithms/name.h b/external/algorithms/core/include/algorithms/name.h new file mode 100644 index 0000000000000000000000000000000000000000..d1668dcf1f82dfc3b039396125158992af1de313 --- /dev/null +++ b/external/algorithms/core/include/algorithms/name.h @@ -0,0 +1,16 @@ +#pragma once + +namespace algorithms { + +// Simple name Mixin providing consistent name API +class NameMixin { +public: + NameMixin(std::string_view name) : m_name{name} {} + std::string_view name() const { return m_name; } + +private: + const std::string m_name; +}; + +} // namespace algorithms + diff --git a/external/algorithms/core/include/algorithms/property.h b/external/algorithms/core/include/algorithms/property.h new file mode 100644 index 0000000000000000000000000000000000000000..e9b3437cf414392290217fcadd74287b0e36c2f6 --- /dev/null +++ b/external/algorithms/core/include/algorithms/property.h @@ -0,0 +1,169 @@ +#pragma once + +#include +#include +#include +#include + +#include + +#include +#include + +namespace algorithms { + +class PropertyError : public Error { +public: + PropertyError(std::string_view msg) : Error{msg, "algorithms::PropertyError"} {} +}; + +// Configuration/property handling +class Configurable { +public: + class PropertyBase; + using PropertyMap = std::map; + + template void setProperty(std::string_view name, T&& value) { + m_props.at(name).set(detail::upcast_type_t(std::forward(value))); + } + template T getProperty(std::string_view name) const { + return static_cast(std::any_cast>(m_props.at(name).get())); + } + const PropertyMap& getProperties() const { return m_props; } + bool hasProperty(std::string_view name) const { + return m_props.count(name) && m_props.at(name).hasValue(); + } + +private: + void registerProperty(PropertyBase& prop) { + if (m_props.count(prop.name())) { + throw PropertyError(fmt::format("Duplicate property name: {}", prop.name())); + } + m_props.emplace(prop.name(), prop); + } + + PropertyMap m_props; + +public: + class PropertyBase { + public: + PropertyBase(std::string_view name) : m_name{name} {} + virtual void set(std::any v) = 0; + virtual std::any get() const = 0; + bool hasValue() const { return m_has_value; } + std::string_view name() const { return m_name; } + + protected: + const std::string m_name; + bool m_has_value = false; + }; + + // A property type that auto-registers itself with the property handler + // Essentially a simplified and const-like version of Gaudi::Property + template class Property : public PropertyBase { + public: + using ValueType = T; + + Property(Configurable* owner, std::string_view name) : PropertyBase{name} { + if (owner) { + owner->registerProperty(*this); + } else { + throw PropertyError( + fmt::format("Attempting to create Property '{}' without valid owner", name)); + } + } + Property(Configurable* owner, std::string_view name, const ValueType& v) + : Property(owner, name) { + set(detail::upcast_type_t(v)); + } + + Property() = delete; + Property(const Property&) = default; + Property& operator=(const Property&) = default; + + // Only settable by explicitly calling the ::set() member functio n + // as we want the Property to mostly act as if it is constant + virtual void set(std::any v) { + m_value = static_cast(std::any_cast>(v)); + m_has_value = true; + } + // virtual getter for use from PropertyBase - use ::value() instead for a quick + // version that does not go through std::any + // Get const reference to the value + virtual std::any get() const { return m_value; } + + // Direct access to the value. Use this one whenever possible (or go through the + // automatic casting) + const ValueType& value() const { return m_value; } + + // automatically cast to T + operator T() const { return m_value; } + + // act as if this is a const T + template bool operator==(const U& rhs) const { return m_value == rhs; } + template bool operator!=(const U& rhs) const { return m_value != rhs; } + template bool operator>(const U& rhs) const { return m_value > rhs; } + template bool operator>=(const U& rhs) const { return m_value >= rhs; } + template bool operator<(const U& rhs) const { return m_value < rhs; } + template bool operator<=(const U& rhs) const { return m_value <= rhs; } + template decltype(auto) operator+(const U& rhs) const { return m_value + rhs; } + template decltype(auto) operator-(const U& rhs) const { return m_value - rhs; } + template decltype(auto) operator*(const U& rhs) const { return m_value * rhs; } + template decltype(auto) operator/(const U& rhs) const { return m_value / rhs; } + + // stl collection helpers if needed + // forced to be templated so we only declare them when used + template decltype(auto) size() const { return value().size(); } + template decltype(auto) length() const { return value().length(); } + template decltype(auto) empty() const { return value().empty(); } + template decltype(auto) clear() { value().clear(); } + template decltype(auto) begin() const { return value().begin(); } + template decltype(auto) end() const { return value().end(); } + template decltype(auto) begin() { return value().begin(); } + template decltype(auto) end() { return value().end(); } + template decltype(auto) operator[](const Arg& arg) const { return value()[arg]; } + template decltype(auto) operator[](const Arg& arg) { return value()[arg]; } + template + decltype(auto) find(const typename U::key_type& key) const { + return value().find(key); + } + template decltype(auto) find(const typename U::key_type& key) { + return value().find(key); + } + template decltype(auto) erase(const Arg& arg) { return value().erase(arg); } + + // In case our property has operator (), delegate the operator() + template + decltype(std::declval()(std::declval()...)) + operator()(Args&&... args) const { + return m_value()(std::forward(args)...); + } + + private: + T m_value; + }; +}; + +// Property mixin, provides all the configuration functionality for +// our algorithms and services +// Currently an alias to Configurable +using PropertyMixin = Configurable; + +} // namespace algorithms + +// operator== overload not needed in C++20 as it will call the member version +#if __cpp_impl_three_way_comparison < 201711 +template +bool operator==(const U& rhs, const algorithms::Configurable::Property& p) { + return p == rhs; +} +#endif +template +std::ostream& operator<<(std::ostream& os, const algorithms::Configurable::Property& p) { + return os << p.value(); +} + +// Make Property formateble +template +struct fmt::formatter> : fmt::formatter {}; +// others needed??? TODO diff --git a/external/algorithms/core/include/algorithms/random.h b/external/algorithms/core/include/algorithms/random.h new file mode 100644 index 0000000000000000000000000000000000000000..09ed60417b3f51211ce84633c7a348ec3df02d48 --- /dev/null +++ b/external/algorithms/core/include/algorithms/random.h @@ -0,0 +1,27 @@ +#pragma once + +#include +#include +#include + +#include +#include + +namespace algorithms { + +class RandomSvc : public LoggedService { +public: + // Initialize the random service + void init(); + + double uniform(double a = 0.0, double b = 1.0) const; + double normal(double mu = 0.0, double sigma = 1.0) const; + +private: + // Standard mersenne_twister_engine seeded with rd() + mutable std::mt19937 rng{std::random_device{}()}; + + ALGORITHMS_DEFINE_LOGGED_SERVICE(RandomSvc) +}; + +} // namespace algorithms diff --git a/external/algorithms/core/include/algorithms/service.h b/external/algorithms/core/include/algorithms/service.h new file mode 100644 index 0000000000000000000000000000000000000000..7fa705ccf94c0614e31803d3a059af4b19d4e417 --- /dev/null +++ b/external/algorithms/core/include/algorithms/service.h @@ -0,0 +1,72 @@ +#pragma once + +#include +#include + +#include +#include + +// Add boilerplate to service class definitions +// - singleton --> no public constructor, no copy constructor, no assigmnent operator +// - constructor should be protected so we can actually inherit from this class if needed +// (mostly needed for the service base class) +#define ALGORITHMS_DEFINE_SERVICE(className) \ +protected: \ + className() : Service(#className) {} \ + \ +public: \ + friend class Service; \ + className(const className&) = delete; \ + void operator=(const className&) = delete; \ + constexpr static const char* kName = #className; + +namespace algorithms { + +// Service service --> keeps track of all services :) +// --> exposes the underlying Configurable object of the service so +// we can configure the services by name in the framework +// boundary plugin +// --> Special service (does not use the Service base class to avoid +// circularity) +class ServiceSvc { +public: + using ServiceMapType = std::map; + static ServiceSvc& instance() { + // This is guaranteed to be thread-safe from C++11 onwards. + static ServiceSvc svc; + return svc; + } + void add(std::string_view name, Configurable* svc) { m_services[name] = svc; } + Configurable* service(std::string_view name) const { return m_services.at(name); } + const ServiceMapType& services() const { return m_services; } + bool active(std::string_view name) const { return m_services.count(name); } + +private: + ServiceSvc() = default; + +public: + ServiceSvc(const ServiceSvc&) = delete; + void operator=(const ServiceSvc&) = delete; + +private: + ServiceMapType m_services; +}; + +// Thread-safe lazy-evaluated minimal service system +// CRTP base class to add the instance method +// This could have been part of DEFINE_SERVICE macro, but I think it is better +// to keep the macro magic to a minimum to maximize transparency +template class Service : public PropertyMixin, public NameMixin { +public: + static SvcType& instance() { + // This is guaranteed to be thread-safe from C++11 onwards. + static SvcType svc; + return svc; + } + // constructor for the service base class registers the service, except + // for the ServiceSvc which is its own thing (avoid circularity) + Service(std::string_view name) : NameMixin{name} { ServiceSvc::instance().add(name, this); } +}; + +} // namespace algorithms + diff --git a/external/algorithms/core/include/algorithms/type_traits.h b/external/algorithms/core/include/algorithms/type_traits.h new file mode 100644 index 0000000000000000000000000000000000000000..db6f99df326d572d6753aedb1e9a2a9dc7e6c97f --- /dev/null +++ b/external/algorithms/core/include/algorithms/type_traits.h @@ -0,0 +1,46 @@ +#pragma once + +#include +#include +#include +#include + +// Useful type traits for generically translation framework specific algorithms into +// using `algorithms' + +namespace algorithms { +// Make it easy to handle the two special data types for algorithms: std::optional +// and std::vector where needed +template struct is_vector : std::false_type {}; +template struct is_vector> : std::true_type {}; +template constexpr bool is_vector_v = is_vector::value; + +template struct is_optional : std::false_type {}; +template struct is_optional> : std::true_type {}; +template constexpr bool is_optional_v = is_optional::value; + +// Get the underlying type for each of the 3 supported cases +template struct data_type { using type = T; }; +template struct data_type> { using type = T; }; +template struct data_type> { using type = T; }; +template using data_type_t = typename data_type::type; + +// Deduce inptu and output value types +template struct deduce_type { + using input_type = gsl::not_null; + using output_type = gsl::not_null; +}; +template struct deduce_type> { + using input_type = const std::vector>; + using output_type = const std::vector>; +}; +template struct deduce_type> { + using input_type = const T*; + using output_type = T*; +}; + +template using input_type_t = typename deduce_type::input_type; +template using output_type_t = typename deduce_type::output_type; + +} // namespace algorithms + diff --git a/external/algorithms/core/src/dummy.cpp b/external/algorithms/core/src/dummy.cpp new file mode 100644 index 0000000000000000000000000000000000000000..5d0bfe6095d4ab1b6e37547e088c21460aa95489 --- /dev/null +++ b/external/algorithms/core/src/dummy.cpp @@ -0,0 +1,86 @@ +#include + +#include + +#include + +using namespace algorithms; + +class testLogger : public LoggerMixin { +public: + testLogger() : LoggerMixin("testLogger") { + std::cout << "should display an info message" << std::endl; + info() << "1" << endmsg; + std::cout << "next debug message should not display" << std::endl; + debug() << "2" << endmsg; + level(LogLevel::kTrace); + std::cout << "Changed the log level to trace, so the following message should display" + << std::endl; + debug() << "3" << endmsg; + std::cout << "error message should display" << std::endl; + error() << "4" << endmsg; + + std::cout << std::endl; + info() << "Checking the existing services:" << endmsg; + for (const auto [key, value] : ServiceSvc::instance().services()) { + info() << " - " << key << endmsg; + } + } +}; + +class propTest : public PropertyMixin, LoggerMixin { +public: + propTest() : LoggerMixin("propTest") {} + void print() const { + info() << "integer property: " << m_int << endmsg; + info() << "foo property: " << m_foo << endmsg; + if (hasProperty("bar")) { + info() << "bar property: " << m_bar << endmsg; + } + if (hasProperty("double")) { + info() << "double (non-def.) property: " << m_double << endmsg; + } + } + +private: + Property m_int{this, "integer", 3}; + Property m_foo{this, "foo", "foo_property_value"}; + // and one without a default + Property m_bar{this, "bar"}; + Property m_double{this, "double"}; +}; + +int dummy() { + testLogger tl; + + std::cout << "test of property, will print set properties. Should only print 2 to start (no bar, " + "as it does not have a default value\n"; + propTest tp; + tp.print(); + tp.setProperty("integer", 10); + tp.print(); + try { + tp.setProperty("foo", 1.); + std::cout << "Should not have been reached, as foo does not exist" << std::endl; + return 1; + } catch (...) { + std::cout << "setting a property that didn't exist failed popperly" << std::endl; + } + std::cout + << "Let's now also set a value for the third and fourth property so it gets printed as well" + << std::endl; + tp.setProperty("bar", "bar_value"); + tp.setProperty("double", 3.1415f); + tp.print(); + + Algorithm, Output>> a{ + "myAlgo", {"int", "double"}, {"moredouble", "variable"}}; + fmt::print("Algo input: {}\n", a.inputNames()); + fmt::print("Algo output: {}\n", a.outputNames()); + // The following won't compile as the number of strings (variable names) and + // input/output tuples don't match + // Algorithm, Output>> a{ + // "myAlgo", {"int", "double"}, {"moredouble", "variable", "bar"}}; + + return 0; +} diff --git a/external/algorithms/core/src/geo.cpp b/external/algorithms/core/src/geo.cpp new file mode 100644 index 0000000000000000000000000000000000000000..ca5e2b68b2d62fb40994e377091f48a92ec6458e --- /dev/null +++ b/external/algorithms/core/src/geo.cpp @@ -0,0 +1,27 @@ +#include + +namespace algorithms { + +void GeoSvc::init(dd4hep::Detector* det) { + if (det) { + info() << "Initializing geometry service from pre-initialized detector" << endmsg; + m_detector = det; + // no detector given, need to self-initialize + } else { + info() << "No external detector provided, self-initializing" << endmsg; + m_detector = &(dd4hep::Detector::getInstance()); + if (m_xml_list.empty()) { + // TODO handle error + } + for (std::string_view name : m_xml_list) { + info() << fmt::format("Loading compact file: {}", "name") << endmsg; + m_detector->fromCompact(std::string(name)); + } + m_detector->volumeManager(); + m_detector->apply("DD4hepVolumeManager", 0, nullptr); + } + // always: instantiate cellIDConverter + m_converter = std::make_unique(*m_detector); +} +} // namespace algorithms + diff --git a/external/algorithms/core/src/random.cpp b/external/algorithms/core/src/random.cpp new file mode 100644 index 0000000000000000000000000000000000000000..16f62934c48d9d71c173486fbf32cfec5f9a651d --- /dev/null +++ b/external/algorithms/core/src/random.cpp @@ -0,0 +1,42 @@ +#include + +#include +#include +#include +#include + +namespace algorithms { + +void RandomSvc::init() { } + +double RandomSvc::uniform(double a, double b) const { + // initialize the random uniform number generator (runif) in a range 0 to 1 + static std::uniform_real_distribution<> runif(a, b); + return runif(rng); +} + +double RandomSvc::normal(double mu, double sigma) const { + // https://en.wikipedia.org/wiki/Box%E2%80%93Muller_transform#Implementation + constexpr double epsilon = std::numeric_limits::epsilon(); + constexpr double two_pi = 2.0 * M_PI; + + static std::uniform_real_distribution<> runif(0.0, 1.0); + + double u1, u2; + do + { + u1 = runif(rng); + } + while (u1 <= epsilon); + u2 = runif(rng); + + //compute z0 and z1 + auto mag = sigma * sqrt(-2.0 * log(u1)); + auto z0 = mag * cos(two_pi * u2) + mu; + //auto z1 = mag * sin(two_pi * u2) + mu; + + return z0; +} + +} // namespace algorithms +