diff --git a/.clang-format b/.clang-format index 5737aca..5c443be 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/.github/workflows/linux-eic-shell.yml b/.github/workflows/linux-eic-shell.yml index 8863a10..54bf320 100644 --- a/.github/workflows/linux-eic-shell.yml +++ b/.github/workflows/linux-eic-shell.yml @@ -18,9 +18,14 @@ jobs: cmake --build build -- install - uses: actions/upload-artifact@v3 with: - name: build-eic-shell + name: install-eic-shell path: install/ if-no-files-found: error + - uses: actions/upload-artifact@v3 + with: + name: build-eic-shell + path: build/ + if-no-files-found: error clang-tidy: runs-on: ubuntu-latest @@ -31,12 +36,12 @@ jobs: - uses: actions/download-artifact@v3 with: name: build-eic-shell - path: install/ + path: build/ - uses: eic/run-cvmfs-osg-eic-shell@main with: platform-release: "jug_xl:nightly" run: | - run-clang-tidy-13 -p build -export-fixes clang_tidy_fixes.yml -extra-arg='-std=c++17' + run-clang-tidy -p build -export-fixes clang_tidy_fixes.yml -extra-arg='-std=c++17' - uses: actions/upload-artifact@v3 with: name: clang-tidy-fixes.yml diff --git a/.gitlab-ci.yml b/.gitlab-ci.yml index 977165b..a1553be 100644 --- a/.gitlab-ci.yml +++ b/.gitlab-ci.yml @@ -12,7 +12,8 @@ variables: ## Other versions to be included in the container. In general these should be "master" ## to keep the container functionally identical to jug_xl:nightly JUGGLER_NPDET_VERSION: "master" - JUGGLER_EICD_VERSION: "master" + JUGGLER_EDM4EIC_VERSION: "main" + JUGGLER_EDM4EIC_REPOSITORYURL: "https://github.com/eic/EDM4eic.git" ## We have: ## - Juggler triggers eic_container on a master pipeline @@ -75,13 +76,13 @@ juggler:local: - CMAKE_CXX_STANDARD: - 17 script: - ## first install EICD to ensure the latest version, then build juggler + ## first install EDM4EIC to ensure the latest version, then build juggler - | - git clone https://eicweb.phy.anl.gov/eic/eicd.git /tmp/eicd - cd /tmp/eicd && git checkout $JUGGLER_EICD_VERSION && cd - - cmake -B /tmp/build -S /tmp/eicd -DCMAKE_CXX_STANDARD=${CMAKE_CXX_STANDARD} -DCMAKE_INSTALL_PREFIX=/usr/local + git clone ${JUGGLER_EDM4EIC_REPOSITORYURL} /tmp/EDM4eic + cd /tmp/EDM4eic && git checkout ${JUGGLER_EDM4EIC_VERSION} && cd - + cmake -B /tmp/build -S /tmp/EDM4eic -DCMAKE_CXX_STANDARD=${CMAKE_CXX_STANDARD} -DCMAKE_INSTALL_PREFIX=/usr/local cmake --build /tmp/build -j40 -- install - rm -rf /tmp/build /tmp/eicd + rm -rf /tmp/build /tmp/EDM4eic - | cmake -Bbuild -S. -DCMAKE_CXX_STANDARD=${CMAKE_CXX_STANDARD} -DCMAKE_INSTALL_PREFIX=/usr/local cmake --build build -j20 @@ -97,7 +98,7 @@ analysis:clang-tidy: - juggler:local script: - | - run-clang-tidy-13 -p build -j20 -export-fixes clang_tidy_fixes.yml -extra-arg='-std=c++17' + run-clang-tidy -p build -j20 -export-fixes clang_tidy_fixes.yml -extra-arg='-std=c++17' artifacts: expire_in: 1 week paths: @@ -130,7 +131,7 @@ version: juggler:default: stage: docker tags: - - docker + - docker-new rules: - if: '$CI_SERVER_HOST == "eicweb.phy.anl.gov"' needs: @@ -145,7 +146,8 @@ juggler:default: --build-arg INTERNAL_TAG=${EIC_DEV_TAG} --build-arg JUGGLER_VERSION=${CI_COMMIT_REF_NAME} --build-arg NPDET_VERSION=${JUGGLER_NPDET_VERSION} - --build-arg EICD_VERSION=${JUGGLER_EICD_VERSION} + --build-arg EDM4EIC_REPOSITORYURL=${JUGGLER_EDM4EIC_REPOSITORYURL} + --build-arg EDM4EIC_VERSION=${JUGGLER_EDM4EIC_VERSION} --build-arg DETECTOR_VERSION=${JUGGLER_DETECTOR_VERSION} --build-arg IP6_VERSION=${JUGGLER_IP6_VERSION} --build-arg JUG_VERSION=juggler-${INTERNAL_TAG}-$(date +%Y-%m-%d_%H-%M-%S)-$(git rev-parse HEAD) diff --git a/CMakeLists.txt b/CMakeLists.txt index 7f17890..d8d743f 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.") @@ -51,7 +51,7 @@ endif() find_package(Microsoft.GSL CONFIG) -find_package(EICD REQUIRED) +find_package(EDM4EIC REQUIRED) find_package(EDM4HEP 0.4.1 REQUIRED) find_package(podio 0.14.1 REQUIRED) @@ -73,11 +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}") -find_library(genfit2 genfit2 /usr/local/lib REQUIRED) -find_path(genfit2_INCLUDE_DIR NAMES GFGbl.h PATHS /usr/local/include ${genfit2}/../include REQUIRED) +## algorithms dependency +add_subdirectory(external/algorithms) find_package(Gaudi) add_subdirectory(JugBase) +add_subdirectory(JugAlgo) add_subdirectory(JugDigi) add_subdirectory(JugFast) add_subdirectory(JugPID) @@ -85,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 0000000..d97bace --- /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 0000000..937b372 --- /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 0000000..bc0c89b --- /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 0000000..77b77c3 --- /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 0000000..2a48f8f --- /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 0000000..4af7933 --- /dev/null +++ b/JugAlgo/src/components/AlgoServiceSvc.cpp @@ -0,0 +1,95 @@ +// 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; + // Always initialize the LogSvc first to ensure proper logging for the others + { + auto& logger = algorithms::LogSvc::instance(); + 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); + } + // loop over all remaining 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) { + ; // Logsvc already initialized, do nothing + } 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()); + } else { + fatal() << "Unknown service encountered, please implement the necessary framework hooks" + << endmsg; + return StatusCode::FAILURE; + } + } + + 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 0000000..959f352 --- /dev/null +++ b/JugAlgo/src/dummy.cpp @@ -0,0 +1,6 @@ +#include +//#include + +namespace { +constexpr int doNothing() { return 1; } +} // namespace diff --git a/JugBase/CMakeLists.txt b/JugBase/CMakeLists.txt index 97bdafb..354bfea 100644 --- a/JugBase/CMakeLists.txt +++ b/JugBase/CMakeLists.txt @@ -19,13 +19,11 @@ gaudi_add_library(JugBase ROOT::Core ROOT::RIO ROOT::Tree DD4hep::DDG4IO DD4hep::DDRec ActsCore ActsPluginDD4hep - ${genfit2} ) target_include_directories(JugBase PUBLIC $ $ - ${genfit2_INCLUDE_DIR} ) target_compile_options(JugBase PRIVATE -Wno-suggest-override) @@ -41,14 +39,12 @@ gaudi_add_module(JugBasePlugins EDM4HEP::edm4hep DD4hep::DDRec ActsCore ActsPluginDD4hep ActsPluginJson - EICD::eicd - ${genfit2} + EDM4EIC::edm4eic ) target_include_directories(JugBasePlugins PUBLIC $ $ - ${genfit2_INCLUDE_DIR} ) target_compile_options(JugBasePlugins PRIVATE -Wno-suggest-override) diff --git a/JugBase/JugBase/IGeoSvc.h b/JugBase/JugBase/IGeoSvc.h index 16d4f18..4ee60fc 100644 --- a/JugBase/JugBase/IGeoSvc.h +++ b/JugBase/JugBase/IGeoSvc.h @@ -22,10 +22,6 @@ namespace Acts { class MagneticFieldProvider; } -namespace genfit { - class DetPlane; -} - /** Geometry service interface. * * \ingroup base @@ -49,13 +45,6 @@ class GAUDI_API IGeoSvc : virtual public IService { virtual const VolumeSurfaceMap& surfaceMap() const = 0; - // Note this hsould return a const& but is just copied for the moment to get around genfit's api - /// Genfit DetPlane map - virtual std::map> getDetPlaneMap() const = 0; - virtual std::map< int64_t, dd4hep::rec::Surface* > getDD4hepSurfaceMap() const =0; - - //virtual std::map< int64_t, dd4hep::rec::Surface* > getDetPlaneMap() const = 0 ; - virtual ~IGeoSvc() {} }; diff --git a/JugBase/JugBase/Utilities/Beam.h b/JugBase/JugBase/Utilities/Beam.h index 9895f56..ebb9438 100644 --- a/JugBase/JugBase/Utilities/Beam.h +++ b/JugBase/JugBase/Utilities/Beam.h @@ -7,7 +7,7 @@ using ROOT::Math::PxPyPzEVector; #include "edm4hep/MCParticleCollection.h" -#include "eicd/ReconstructedParticleCollection.h" +#include "edm4eic/ReconstructedParticleCollection.h" namespace Jug::Base::Beam { @@ -53,7 +53,7 @@ namespace Jug::Base::Beam { return find_first_with_status_pdg(mcparts, {1}, {11}); } - inline auto find_first_scattered_electron(const eicd::ReconstructedParticleCollection& rcparts) { + inline auto find_first_scattered_electron(const edm4eic::ReconstructedParticleCollection& rcparts) { return find_first_with_pdg(rcparts, {11}); } diff --git a/JugBase/src/components/GeoSvc.cpp b/JugBase/src/components/GeoSvc.cpp index 2eddc4b..c6d3554 100644 --- a/JugBase/src/components/GeoSvc.cpp +++ b/JugBase/src/components/GeoSvc.cpp @@ -18,20 +18,6 @@ #include "Acts/MagneticField/MagneticFieldContext.hpp" #include "Acts/Surfaces/PlaneSurface.hpp" -// genfit -#include "ConstField.h" -#include "DAF.h" -#include "Exception.h" -#include "FieldManager.h" -#include "KalmanFitterRefTrack.h" -#include "StateOnPlane.h" -#include "Track.h" -#include "TrackPoint.h" -#include "MaterialEffects.h" -#include "RKTrackRep.h" -#include "TGeoMaterialInterface.h" -#include "PlanarMeasurement.h" - static const std::map s_msgMap = { {MSG::DEBUG, Acts::Logging::DEBUG}, {MSG::VERBOSE, Acts::Logging::VERBOSE}, @@ -113,29 +99,6 @@ StatusCode GeoSvc::initialize() { m_log << MSG::INFO << "DD4Hep geometry SUCCESSFULLY built" << endmsg; } - // Genfit - genfit::FieldManager::getInstance()->init(new genfit::ConstField( - 0., 0., this->centralMagneticField() * 10.0)); // gentfit uses kilo-Gauss - genfit::MaterialEffects::getInstance()->init(new genfit::TGeoMaterialInterface()); - - // create a list of all surfaces in the detector: - dd4hep::rec::SurfaceManager surfMan( *m_dd4hepGeo ) ; - debug() << " surface manager " << endmsg; - const auto* const sM = surfMan.map("tracker") ; - if (sM != nullptr) { - debug() << " surface map size: " << sM->size() << endmsg; - // setup dd4hep surface map - //for( dd4hep::rec::SurfaceMap::const_iterator it = sM->begin() ; it != sM->end() ; ++it ){ - for( const auto& [id, s] : *sM) { - //dd4hep::rec::Surface* surf = s ; - m_surfaceMap[ id ] = dynamic_cast(s) ; - debug() << " surface : " << *s << endmsg; - m_detPlaneMap[id] = std::shared_ptr( - new genfit::DetPlane({s->origin().x(), s->origin().y(), s->origin().z()}, {s->u().x(), s->u().y(), s->u().z()}, - {s->v().x(), s->v().y(), s->v().z()})); - } - } - // Set ACTS logging level auto im = s_msgMap.find(msgLevel()); if (im != s_msgMap.end()) { diff --git a/JugBase/src/components/GeoSvc.h b/JugBase/src/components/GeoSvc.h index d8ea303..d992de0 100644 --- a/JugBase/src/components/GeoSvc.h +++ b/JugBase/src/components/GeoSvc.h @@ -57,12 +57,6 @@ class GeoSvc : public extends { */ dd4hep::Detector* m_dd4hepGeo = nullptr; - /// DD4hep surface map - std::map< int64_t, dd4hep::rec::Surface* > m_surfaceMap ; - - /// Genfit DetPlane map - std::map< int64_t, std::shared_ptr > m_detPlaneMap ; - /// ACTS Logging Level Acts::Logging::Level m_actsLoggingLevel = Acts::Logging::INFO; @@ -144,11 +138,6 @@ class GeoSvc : public extends { } virtual const VolumeSurfaceMap& surfaceMap() const { return m_surfaces; } - - // Note this hsould return a const& but is just copied for the moment to get around genfit's api - virtual std::map> getDetPlaneMap() const { return m_detPlaneMap; } - - virtual std::map< int64_t, dd4hep::rec::Surface* > getDD4hepSurfaceMap() const { return m_surfaceMap ;} }; inline std::shared_ptr GeoSvc::trackingGeometry() const diff --git a/JugBase/src/components/PodioOutput.cpp b/JugBase/src/components/PodioOutput.cpp index 2793878..fe2f97e 100644 --- a/JugBase/src/components/PodioOutput.cpp +++ b/JugBase/src/components/PodioOutput.cpp @@ -76,10 +76,6 @@ void PodioOutput::resetBranches(const std::vector>& collections) { - // collectionID, collection type, subset collection - auto* collectionInfo = new std::vector>(); - collectionInfo->reserve(collections.size()); - for (const auto& [collName, collBuffers] : collections) { auto buffers = collBuffers->getBuffers(); auto* data = buffers.data; @@ -116,22 +112,22 @@ void PodioOutput::createBranches(const std::vectorgetCollectionIDs()->collectionID(collName); const auto collType = collBuffers->getValueTypeName() + "Collection"; - collectionInfo->emplace_back(collID, std::move(collType), collBuffers->isSubsetCollection()); + m_collectionInfo.emplace_back(collID, std::move(collType), collBuffers->isSubsetCollection()); debug() << isOn << " Registering collection " << collClassName << " " << collName.c_str() << " containing type " << className << endmsg; collBuffers->prepareForWrite(); } - - m_metadatatree->Branch("CollectionTypeInfo", collectionInfo); } StatusCode PodioOutput::execute() { // for now assume identical content for every event // register for writing if (m_firstEvent) { + m_collectionInfo.clear(); createBranches(m_podioDataSvc->getCollections()); createBranches(m_podioDataSvc->getReadCollections()); + m_metadatatree->Branch("CollectionTypeInfo", &m_collectionInfo); } else { resetBranches(m_podioDataSvc->getCollections()); resetBranches(m_podioDataSvc->getReadCollections()); diff --git a/JugBase/src/components/PodioOutput.h b/JugBase/src/components/PodioOutput.h index 5bb6408..2dba4dd 100644 --- a/JugBase/src/components/PodioOutput.h +++ b/JugBase/src/components/PodioOutput.h @@ -58,6 +58,7 @@ class PodioOutput : public GaudiAlgorithm { gsl::owner m_colMDtree; /// The stored collections std::vector m_storedCollections; + std::vector> m_collectionInfo; }; #endif diff --git a/JugDigi/CMakeLists.txt b/JugDigi/CMakeLists.txt index 64a0749..41d082b 100644 --- a/JugDigi/CMakeLists.txt +++ b/JugDigi/CMakeLists.txt @@ -14,7 +14,7 @@ gaudi_add_module(JugDigiPlugins JugBase ROOT::Core ROOT::RIO ROOT::Tree EDM4HEP::edm4hep - EICD::eicd + EDM4EIC::edm4eic ) target_include_directories(JugDigiPlugins PUBLIC diff --git a/JugDigi/src/components/CalorimeterBirksCorr.cpp b/JugDigi/src/components/CalorimeterBirksCorr.cpp index b42c1dd..392f1bf 100644 --- a/JugDigi/src/components/CalorimeterBirksCorr.cpp +++ b/JugDigi/src/components/CalorimeterBirksCorr.cpp @@ -21,8 +21,8 @@ // Event Model related classes #include "edm4hep/SimCalorimeterHitCollection.h" -#include "eicd/RawCalorimeterHitCollection.h" -#include "eicd/RawCalorimeterHitData.h" +#include "edm4eic/RawCalorimeterHitCollection.h" +#include "edm4eic/RawCalorimeterHitData.h" using namespace Gaudi::Units; diff --git a/JugDigi/src/components/CalorimeterHitDigi.cpp b/JugDigi/src/components/CalorimeterHitDigi.cpp index 6224189..a108d76 100644 --- a/JugDigi/src/components/CalorimeterHitDigi.cpp +++ b/JugDigi/src/components/CalorimeterHitDigi.cpp @@ -31,8 +31,8 @@ // Event Model related classes #include "edm4hep/SimCalorimeterHitCollection.h" -#include "eicd/RawCalorimeterHitCollection.h" -#include "eicd/RawCalorimeterHitData.h" +#include "edm4eic/RawCalorimeterHitCollection.h" +#include "edm4eic/RawCalorimeterHitData.h" using namespace Gaudi::Units; @@ -49,6 +49,8 @@ namespace Jug::Digi { // 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}; @@ -80,7 +82,7 @@ namespace Jug::Digi { DataHandle m_inputHitCollection{ "inputHitCollection", Gaudi::DataHandle::Reader, this}; - DataHandle m_outputHitCollection{ + DataHandle m_outputHitCollection{ "outputHitCollection", Gaudi::DataHandle::Writer, this}; // ill-formed: using GaudiAlgorithm::GaudiAlgorithm; @@ -173,13 +175,16 @@ namespace Jug::Digi { const double eDep = ahit.getEnergy(); // apply additional calorimeter noise to corrected energy deposit - const double eResRel = (eDep > 1e-6) - ? 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 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 + m_corrMeanScale * eDep * (1. + eResRel) / dyRangeADC * m_capADC); + 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()) { @@ -189,7 +194,7 @@ namespace Jug::Digi { } const long long tdc = std::llround((time + m_normDist() * tRes) * stepTDC); - eicd::RawCalorimeterHit rawhit( + edm4eic::RawCalorimeterHit rawhit( ahit.getCellID(), (adc > m_capADC.value() ? m_capADC.value() : adc), tdc @@ -233,18 +238,18 @@ namespace Jug::Digi { } } - double eResRel = 0.; // safety check - if (edep > 1e-6) { - eResRel = m_normDist() * eRes[0] / std::sqrt(edep) + - m_normDist() * eRes[1] + - m_normDist() * eRes[2] / edep; - } + 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); - eicd::RawCalorimeterHit rawhit( + edm4eic::RawCalorimeterHit rawhit( id, (adc > m_capADC.value() ? m_capADC.value() : adc), tdc diff --git a/JugDigi/src/components/PhotoMultiplierDigi.cpp b/JugDigi/src/components/PhotoMultiplierDigi.cpp index c5a5480..2fde77d 100644 --- a/JugDigi/src/components/PhotoMultiplierDigi.cpp +++ b/JugDigi/src/components/PhotoMultiplierDigi.cpp @@ -23,7 +23,7 @@ #include "JugBase/DataHandle.h" // Event Model related classes -#include "eicd/RawPMTHitCollection.h" +#include "edm4eic/RawPMTHitCollection.h" #include "edm4hep/MCParticleCollection.h" #include "edm4hep/SimTrackerHitCollection.h" @@ -41,7 +41,7 @@ class PhotoMultiplierDigi : public GaudiAlgorithm public: DataHandle m_inputHitCollection{"inputHitCollection", Gaudi::DataHandle::Reader, this}; - DataHandle + DataHandle m_outputHitCollection{"outputHitCollection", Gaudi::DataHandle::Writer, this}; Gaudi::Property>> u_quantumEfficiency{this, "quantumEfficiency", {{2.6*eV, 0.3}, {7.0*eV, 0.3}}}; @@ -124,7 +124,7 @@ class PhotoMultiplierDigi : public GaudiAlgorithm // build hit for (auto &it : hit_groups) { for (auto &data : it.second) { - eicd::RawPMTHit hit{ + edm4eic::RawPMTHit hit{ it.first, static_cast(data.signal), static_cast(data.time/(m_timeStep/ns))}; diff --git a/JugDigi/src/components/SiliconTrackerDigi.cpp b/JugDigi/src/components/SiliconTrackerDigi.cpp index 4789ce2..d37bfec 100644 --- a/JugDigi/src/components/SiliconTrackerDigi.cpp +++ b/JugDigi/src/components/SiliconTrackerDigi.cpp @@ -16,8 +16,8 @@ // edm4hep's tracker hit is the input collectiopn #include "edm4hep/MCParticle.h" #include "edm4hep/SimTrackerHitCollection.h" -// eicd's RawTrackerHit is the output -#include "eicd/RawTrackerHitCollection.h" +// edm4eic's RawTrackerHit is the output +#include "edm4eic/RawTrackerHitCollection.h" namespace Jug::Digi { @@ -32,7 +32,7 @@ class SiliconTrackerDigi : public GaudiAlgorithm { Rndm::Numbers m_gaussDist; DataHandle m_inputHitCollection{"inputHitCollection", Gaudi::DataHandle::Reader, this}; - DataHandle m_outputHitCollection{"outputHitCollection", Gaudi::DataHandle::Writer, + DataHandle m_outputHitCollection{"outputHitCollection", Gaudi::DataHandle::Writer, this}; public: @@ -57,7 +57,7 @@ class SiliconTrackerDigi : public GaudiAlgorithm { const auto* const simhits = m_inputHitCollection.get(); // Create output collections auto* rawhits = m_outputHitCollection.createAndPut(); - // eicd::RawTrackerHitCollection* rawHitCollection = new eicd::RawTrackerHitCollection(); + // edm4eic::RawTrackerHitCollection* rawHitCollection = new edm4eic::RawTrackerHitCollection(); std::map cell_hit_map; for (const auto& ahit : *simhits) { if (msgLevel(MSG::DEBUG)) { @@ -82,7 +82,7 @@ class SiliconTrackerDigi : public GaudiAlgorithm { } if (cell_hit_map.count(ahit.getCellID()) == 0) { cell_hit_map[ahit.getCellID()] = rawhits->size(); - eicd::RawTrackerHit rawhit(ahit.getCellID(), + edm4eic::RawTrackerHit rawhit(ahit.getCellID(), ahit.getMCParticle().getTime() * 1e6 + m_gaussDist() * 1e3, // ns->fs std::llround(ahit.getEDep() * 1e6)); rawhits->push_back(rawhit); diff --git a/JugFast/CMakeLists.txt b/JugFast/CMakeLists.txt index 96d3ccf..91b9c9b 100644 --- a/JugFast/CMakeLists.txt +++ b/JugFast/CMakeLists.txt @@ -14,7 +14,7 @@ gaudi_add_module(JugFastPlugins JugBase ROOT::Core ROOT::RIO ROOT::Tree EDM4HEP::edm4hep - EICD::eicd + EDM4EIC::edm4eic ) target_include_directories(JugFastPlugins PUBLIC diff --git a/JugFast/src/components/ClusterMerger.cpp b/JugFast/src/components/ClusterMerger.cpp index 30df222..45b47b1 100644 --- a/JugFast/src/components/ClusterMerger.cpp +++ b/JugFast/src/components/ClusterMerger.cpp @@ -15,9 +15,9 @@ #include "JugBase/DataHandle.h" // Event Model related classes -#include "eicd/ClusterCollection.h" -#include "eicd/MCRecoClusterParticleAssociationCollection.h" -#include "eicd/vector_utils.h" +#include "edm4eic/ClusterCollection.h" +#include "edm4eic/MCRecoClusterParticleAssociationCollection.h" +#include "edm4eic/vector_utils.h" using namespace Gaudi::Units; @@ -31,11 +31,11 @@ namespace Jug::Fast { class ClusterMerger : public GaudiAlgorithm { private: // Input - DataHandle m_inputClusters{"InputClusters", Gaudi::DataHandle::Reader, this}; - DataHandle m_inputAssociations{"InputAssociations", Gaudi::DataHandle::Reader, this}; + DataHandle m_inputClusters{"InputClusters", Gaudi::DataHandle::Reader, this}; + DataHandle m_inputAssociations{"InputAssociations", Gaudi::DataHandle::Reader, this}; // Output - DataHandle m_outputClusters{"OutputClusters", Gaudi::DataHandle::Writer, this}; - DataHandle m_outputAssociations{"OutputAssociations", Gaudi::DataHandle::Writer, this}; + DataHandle m_outputClusters{"OutputClusters", Gaudi::DataHandle::Writer, this}; + DataHandle m_outputAssociations{"OutputAssociations", Gaudi::DataHandle::Writer, this}; public: ClusterMerger(const std::string& name, ISvcLocator* svcLoc) : GaudiAlgorithm(name, svcLoc) { @@ -106,7 +106,7 @@ class ClusterMerger : public GaudiAlgorithm { float energyError = 0; float time = 0; int nhits = 0; - eicd::Vector3f position; + auto position = new_clus.getPosition(); for (const auto& clus : clusters) { if (msgLevel(MSG::DEBUG)) { debug() << " --> Adding cluster with energy: " << clus.getEnergy() << endmsg; @@ -142,12 +142,12 @@ class ClusterMerger : public GaudiAlgorithm { } // get a map of MCParticle index--> std::vector for clusters that belong together - std::map> indexedClusterLists( - const eicd::ClusterCollection& clusters, - const eicd::MCRecoClusterParticleAssociationCollection& associations + std::map> indexedClusterLists( + const edm4eic::ClusterCollection& clusters, + const edm4eic::MCRecoClusterParticleAssociationCollection& associations ) const { - std::map> matched = {}; + std::map> matched = {}; // loop over clusters for (const auto& cluster : clusters) { diff --git a/JugFast/src/components/InclusiveKinematicsTruth.cpp b/JugFast/src/components/InclusiveKinematicsTruth.cpp index 0f35704..48f7e34 100644 --- a/JugFast/src/components/InclusiveKinematicsTruth.cpp +++ b/JugFast/src/components/InclusiveKinematicsTruth.cpp @@ -20,9 +20,9 @@ using ROOT::Math::PxPyPzEVector; // Event Model related classes #include "edm4hep/MCParticleCollection.h" -#include "eicd/InclusiveKinematicsCollection.h" +#include "edm4eic/InclusiveKinematicsCollection.h" -#include "eicd/vector_utils.h" +#include "edm4eic/vector_utils.h" namespace Jug::Fast { @@ -32,7 +32,7 @@ class InclusiveKinematicsTruth : public GaudiAlgorithm { "inputMCParticles", Gaudi::DataHandle::Reader, this}; - DataHandle m_outputInclusiveKinematicsCollection{ + DataHandle m_outputInclusiveKinematicsCollection{ "outputInclusiveKinematics", Gaudi::DataHandle::Writer, this}; @@ -89,7 +89,7 @@ class InclusiveKinematicsTruth : public GaudiAlgorithm { return StatusCode::SUCCESS; } const auto ei_p = ei_coll[0].getMomentum(); - const auto ei_p_mag = eicd::magnitude(ei_p); + const auto ei_p_mag = edm4eic::magnitude(ei_p); const auto ei_mass = m_electron; const PxPyPzEVector ei(ei_p.x, ei_p.y, ei_p.z, std::hypot(ei_p_mag, ei_mass)); @@ -102,7 +102,7 @@ class InclusiveKinematicsTruth : public GaudiAlgorithm { return StatusCode::SUCCESS; } const auto pi_p = pi_coll[0].getMomentum(); - const auto pi_p_mag = eicd::magnitude(pi_p); + const auto pi_p_mag = edm4eic::magnitude(pi_p); const auto pi_mass = pi_coll[0].getPDG() == 2212 ? m_proton : m_neutron; const PxPyPzEVector pi(pi_p.x, pi_p.y, pi_p.z, std::hypot(pi_p_mag, pi_mass)); @@ -119,7 +119,7 @@ class InclusiveKinematicsTruth : public GaudiAlgorithm { return StatusCode::SUCCESS; } const auto ef_p = ef_coll[0].getMomentum(); - const auto ef_p_mag = eicd::magnitude(ef_p); + const auto ef_p_mag = edm4eic::magnitude(ef_p); const auto ef_mass = m_electron; const PxPyPzEVector ef(ef_p.x, ef_p.y, ef_p.z, std::hypot(ef_p_mag, ef_mass)); diff --git a/JugFast/src/components/MC2SmearedParticle.cpp b/JugFast/src/components/MC2SmearedParticle.cpp index b3e691d..6968e37 100644 --- a/JugFast/src/components/MC2SmearedParticle.cpp +++ b/JugFast/src/components/MC2SmearedParticle.cpp @@ -13,14 +13,14 @@ // Event Model related classes #include "edm4hep/MCParticleCollection.h" -#include "eicd/ReconstructedParticleCollection.h" +#include "edm4eic/ReconstructedParticleCollection.h" namespace Jug::Fast { class MC2SmearedParticle : public GaudiAlgorithm { private: DataHandle m_inputMCParticles{"MCParticles", Gaudi::DataHandle::Reader, this}; - DataHandle m_outputParticles{"SmearedReconstructedParticles", + DataHandle m_outputParticles{"SmearedReconstructedParticles", Gaudi::DataHandle::Writer, this}; Rndm::Numbers m_gaussDist; Gaudi::Property m_smearing{this, "smearing", 0.01 /* 1 percent*/}; @@ -62,7 +62,7 @@ class MC2SmearedParticle : public GaudiAlgorithm { const auto pgen = std::hypot(pvec.x, pvec.y, pvec.z); const auto momentum = pgen * m_gaussDist(); // make sure we keep energy consistent - using MomType = decltype(eicd::ReconstructedParticle().getMomentum().x); + using MomType = decltype(edm4eic::ReconstructedParticle().getMomentum().x); const MomType energy = std::sqrt(p.getEnergy() * p.getEnergy() - pgen * pgen + momentum * momentum); const MomType px = p.getMomentum().x * momentum / pgen; const MomType py = p.getMomentum().y * momentum / pgen; diff --git a/JugFast/src/components/MatchClusters.cpp b/JugFast/src/components/MatchClusters.cpp index 547d8af..2a888b3 100644 --- a/JugFast/src/components/MatchClusters.cpp +++ b/JugFast/src/components/MatchClusters.cpp @@ -20,12 +20,12 @@ // Event Model related classes #include "edm4hep/MCParticleCollection.h" -#include "eicd/ClusterCollection.h" -#include "eicd/MCRecoClusterParticleAssociationCollection.h" -#include "eicd/MCRecoParticleAssociationCollection.h" -#include "eicd/ReconstructedParticleCollection.h" -#include "eicd/TrackParametersCollection.h" -#include "eicd/vector_utils.h" +#include "edm4eic/ClusterCollection.h" +#include "edm4eic/MCRecoClusterParticleAssociationCollection.h" +#include "edm4eic/MCRecoParticleAssociationCollection.h" +#include "edm4eic/ReconstructedParticleCollection.h" +#include "edm4eic/TrackParametersCollection.h" +#include "edm4eic/vector_utils.h" namespace Jug::Fast { @@ -33,19 +33,19 @@ class MatchClusters : public GaudiAlgorithm { private: // input data DataHandle m_inputMCParticles{"MCParticles", Gaudi::DataHandle::Reader, this}; - DataHandle m_inputParticles{"ReconstructedChargedParticles", + DataHandle m_inputParticles{"ReconstructedChargedParticles", Gaudi::DataHandle::Reader, this}; - DataHandle m_inputParticlesAssoc{"ReconstructedChargedParticlesAssoc", + DataHandle m_inputParticlesAssoc{"ReconstructedChargedParticlesAssoc", Gaudi::DataHandle::Reader, this}; Gaudi::Property> m_inputClusters{this, "inputClusters", {}, "Clusters to be aggregated"}; Gaudi::Property> m_inputClustersAssoc{this, "inputClustersAssoc", {}, "Cluster associations to be aggregated"}; - std::vector*> m_inputClustersCollections; - std::vector*> m_inputClustersAssocCollections; + std::vector*> m_inputClustersCollections; + std::vector*> m_inputClustersAssocCollections; // output data - DataHandle m_outputParticles{"ReconstructedParticles", + DataHandle m_outputParticles{"ReconstructedParticles", Gaudi::DataHandle::Writer, this}; - DataHandle m_outputParticlesAssoc{"ReconstructedParticlesAssoc", + DataHandle m_outputParticlesAssoc{"ReconstructedParticlesAssoc", Gaudi::DataHandle::Writer, this}; public: @@ -183,32 +183,32 @@ class MatchClusters : public GaudiAlgorithm { } private: - std::vector*> getClusterCollections(const std::vector& cols) { - std::vector*> ret; + std::vector*> getClusterCollections(const std::vector& cols) { + std::vector*> ret; for (const auto& colname : cols) { debug() << "initializing cluster collection: " << colname << endmsg; - ret.push_back(new DataHandle{colname, Gaudi::DataHandle::Reader, this}); + ret.push_back(new DataHandle{colname, Gaudi::DataHandle::Reader, this}); } return ret; } - std::vector*> getClusterAssociations(const std::vector& cols) { - std::vector*> ret; + std::vector*> getClusterAssociations(const std::vector& cols) { + std::vector*> ret; for (const auto& colname : cols) { debug() << "initializing cluster association collection: " << colname << endmsg; - ret.push_back(new DataHandle{colname, Gaudi::DataHandle::Reader, this}); + ret.push_back(new DataHandle{colname, Gaudi::DataHandle::Reader, this}); } return ret; } // get a map of mcID --> cluster // input: cluster_collections --> list of handles to all cluster collections - std::map + std::map indexedClusters( - const std::vector*>& cluster_collections, - const std::vector*>& associations_collections + const std::vector*>& cluster_collections, + const std::vector*>& associations_collections ) const { - std::map matched = {}; + std::map matched = {}; // loop over cluster collections for (const auto& cluster_handle : cluster_collections) { @@ -266,14 +266,14 @@ class MatchClusters : public GaudiAlgorithm { // reconstruct a neutral cluster // (for now assuming the vertex is at (0,0,0)) - eicd::ReconstructedParticle reconstruct_neutral(const eicd::Cluster& clus, const double mass, + edm4eic::ReconstructedParticle reconstruct_neutral(const edm4eic::Cluster& clus, const double mass, const int32_t pdg) const { const float energy = clus.getEnergy(); const float p = energy < mass ? 0 : std::sqrt(energy * energy - mass * mass); const auto position = clus.getPosition(); - const auto momentum = p * (position / eicd::magnitude(position)); + const auto momentum = p * (position / edm4eic::magnitude(position)); // setup our particle - eicd::MutableReconstructedParticle part; + edm4eic::MutableReconstructedParticle part; part.setMomentum(momentum); part.setPDG(pdg); part.setCharge(0); diff --git a/JugFast/src/components/ParticlesWithTruthPID.cpp b/JugFast/src/components/ParticlesWithTruthPID.cpp index c09a18e..40177ae 100644 --- a/JugFast/src/components/ParticlesWithTruthPID.cpp +++ b/JugFast/src/components/ParticlesWithTruthPID.cpp @@ -16,21 +16,21 @@ // Event Model related classes #include "edm4hep/MCParticleCollection.h" -#include "eicd/MCRecoParticleAssociationCollection.h" -#include "eicd/ReconstructedParticleCollection.h" -#include "eicd/TrackParametersCollection.h" -#include "eicd/vector_utils.h" +#include "edm4eic/MCRecoParticleAssociationCollection.h" +#include "edm4eic/ReconstructedParticleCollection.h" +#include "edm4eic/TrackParametersCollection.h" +#include "edm4eic/vector_utils.h" namespace Jug::Fast { class ParticlesWithTruthPID : public GaudiAlgorithm { private: DataHandle m_inputTruthCollection{"inputMCParticles", Gaudi::DataHandle::Reader, this}; - DataHandle m_inputTrackCollection{"inputTrackParameters", Gaudi::DataHandle::Reader, + DataHandle m_inputTrackCollection{"inputTrackParameters", Gaudi::DataHandle::Reader, this}; - DataHandle m_outputParticleCollection{"ReconstructedParticles", + DataHandle m_outputParticleCollection{"ReconstructedParticles", Gaudi::DataHandle::Writer, this}; - DataHandle m_outputAssocCollection{"MCRecoParticleAssociation", + DataHandle m_outputAssocCollection{"MCRecoParticleAssociation", Gaudi::DataHandle::Writer, this}; // Matching momentum tolerance requires 10% by default; @@ -63,7 +63,7 @@ class ParticlesWithTruthPID : public GaudiAlgorithm { const double sinPhiOver2Tolerance = sin(0.5 * m_phiTolerance); std::vector consumed(mc.size(), false); for (const auto& trk : tracks) { - const auto mom = eicd::sphericalToVector(1.0 / std::abs(trk.getQOverP()), trk.getTheta(), trk.getPhi()); + const auto mom = edm4eic::sphericalToVector(1.0 / std::abs(trk.getQOverP()), trk.getTheta(), trk.getPhi()); const auto charge_rec = trk.getCharge(); // utility variables for matching int best_match = -1; @@ -81,11 +81,11 @@ class ParticlesWithTruthPID : public GaudiAlgorithm { const auto p_mag = std::hypot(p.x, p.y, p.z); const auto p_phi = std::atan2(p.y, p.x); const auto p_eta = std::atanh(p.z / p_mag); - const double dp_rel = std::abs((eicd::magnitude(mom) - p_mag) / p_mag); + const double dp_rel = std::abs((edm4eic::magnitude(mom) - p_mag) / p_mag); // check the tolerance for sin(dphi/2) to avoid the hemisphere problem and allow // for phi rollovers - const double dsphi = std::abs(sin(0.5 * (eicd::angleAzimuthal(mom) - p_phi))); - const double deta = std::abs((eicd::eta(mom) - p_eta)); + const double dsphi = std::abs(sin(0.5 * (edm4eic::angleAzimuthal(mom) - p_phi))); + const double deta = std::abs((edm4eic::eta(mom) - p_eta)); if (dp_rel < m_pRelativeTolerance && deta < m_etaTolerance && dsphi < sinPhiOver2Tolerance) { const double delta = @@ -112,7 +112,7 @@ class ParticlesWithTruthPID : public GaudiAlgorithm { mass = mcpart.getMass(); } rec_part.setType(static_cast(best_match >= 0 ? 0 : -1)); // @TODO: determine type codes - rec_part.setEnergy(std::hypot(eicd::magnitude(mom), mass)); + rec_part.setEnergy(std::hypot(edm4eic::magnitude(mom), mass)); rec_part.setMomentum(mom); rec_part.setReferencePoint(referencePoint); rec_part.setCharge(charge_rec); @@ -133,20 +133,20 @@ class ParticlesWithTruthPID : public GaudiAlgorithm { if (best_match > 0) { const auto& mcpart = mc[best_match]; debug() << fmt::format("Matched track with MC particle {}\n", best_match) << endmsg; - debug() << fmt::format(" - Track: (mom: {}, theta: {}, phi: {}, charge: {})", eicd::magnitude(mom), - eicd::anglePolar(mom), eicd::angleAzimuthal(mom), charge_rec) + debug() << fmt::format(" - Track: (mom: {}, theta: {}, phi: {}, charge: {})", edm4eic::magnitude(mom), + edm4eic::anglePolar(mom), edm4eic::angleAzimuthal(mom), charge_rec) << endmsg; const auto& p = mcpart.getMomentum(); - const auto p_mag = eicd::magnitude(p); - const auto p_phi = eicd::angleAzimuthal(p); - const auto p_theta = eicd::anglePolar(p); + const auto p_mag = edm4eic::magnitude(p); + const auto p_phi = edm4eic::angleAzimuthal(p); + const auto p_theta = edm4eic::anglePolar(p); debug() << fmt::format(" - MC particle: (mom: {}, theta: {}, phi: {}, charge: {}, type: {}", p_mag, p_theta, p_phi, mcpart.getCharge(), mcpart.getPDG()) << endmsg; } else { debug() << fmt::format("Did not find a good match for track \n") << endmsg; - debug() << fmt::format(" - Track: (mom: {}, theta: {}, phi: {}, charge: {})", eicd::magnitude(mom), - eicd::anglePolar(mom), eicd::angleAzimuthal(mom), charge_rec) + debug() << fmt::format(" - Track: (mom: {}, theta: {}, phi: {}, charge: {})", edm4eic::magnitude(mom), + edm4eic::anglePolar(mom), edm4eic::angleAzimuthal(mom), charge_rec) << endmsg; } } diff --git a/JugFast/src/components/ScatteredElectronFinder.cpp b/JugFast/src/components/ScatteredElectronFinder.cpp new file mode 100644 index 0000000..8ef2794 --- /dev/null +++ b/JugFast/src/components/ScatteredElectronFinder.cpp @@ -0,0 +1,72 @@ +// SPDX-License-Identifier: LGPL-3.0-or-later +// Copyright (C) 2022 Wouter Deconinck + +#include "Gaudi/Algorithm.h" +#include "GaudiAlg/GaudiTool.h" +#include "GaudiAlg/Producer.h" +#include "GaudiAlg/Transformer.h" +#include "GaudiKernel/RndmGenerators.h" +#include "GaudiKernel/PhysicalConstants.h" +#include +#include + +#include "JugBase/IParticleSvc.h" +#include "JugBase/DataHandle.h" + +#include "JugBase/Utilities/Beam.h" + +#include "Math/Vector4D.h" +using ROOT::Math::PxPyPzEVector; + +// Event Model related classes +#include "edm4hep/MCParticleCollection.h" + +namespace Jug::Fast { + +class ScatteredElectronFinder : public GaudiAlgorithm { +private: + DataHandle m_inputMCParticleCollection{ + "inputMCParticles", + Gaudi::DataHandle::Reader, + this}; + DataHandle m_outputMCScatteredElectron{ + "outputMCScatteredElectron", + Gaudi::DataHandle::Writer, + this}; + +public: + ScatteredElectronFinder(const std::string& name, ISvcLocator* svcLoc) + : GaudiAlgorithm(name, svcLoc) { + declareProperty("inputMCParticles", m_inputMCParticleCollection, "MCParticles"); + declareProperty("outputMCScatteredElectron", m_outputMCScatteredElectron, "MCScatteredElectron"); + } + + StatusCode initialize() override { + return GaudiAlgorithm::initialize(); + } + + StatusCode execute() override { + // input collections + const auto& mcparts = *(m_inputMCParticleCollection.get()); + // output collection + auto& out_electron = *(m_outputMCScatteredElectron.createAndPut()); + out_electron.setSubsetCollection(); + + // Determine scattered electron + // + // Currently taken as first status==1 electron in HEPMC record, + // which seems to be correct based on a cursory glance at the + // Pythia8 output. In the future, it may be better to trace back + // each final-state electron and see which one originates from + // the beam. + const auto ef_coll = Jug::Base::Beam::find_first_scattered_electron(mcparts); + out_electron.push_back(ef_coll.front()); + + return StatusCode::SUCCESS; + } +}; + +// NOLINTNEXTLINE(cppcoreguidelines-avoid-non-const-global-variables) +DECLARE_COMPONENT(ScatteredElectronFinder) + +} // namespace Jug::Fast diff --git a/JugFast/src/components/SmearedFarForwardParticles.cpp b/JugFast/src/components/SmearedFarForwardParticles.cpp index 629c1d6..bdebe27 100644 --- a/JugFast/src/components/SmearedFarForwardParticles.cpp +++ b/JugFast/src/components/SmearedFarForwardParticles.cpp @@ -15,9 +15,9 @@ // Event Model related classes #include "edm4hep/MCParticleCollection.h" -#include "eicd/MCRecoParticleAssociationCollection.h" -#include "eicd/ReconstructedParticleCollection.h" -#include "eicd/vector_utils.h" +#include "edm4eic/MCRecoParticleAssociationCollection.h" +#include "edm4eic/ReconstructedParticleCollection.h" +#include "edm4eic/vector_utils.h" namespace { enum DetectorTags { kTagB0 = 1, kTagRP = 2, kTagOMD = 3, kTagZDC = 4 }; @@ -28,9 +28,9 @@ namespace Jug::Fast { class SmearedFarForwardParticles : public GaudiAlgorithm { private: DataHandle m_inputMCParticles{"inputMCParticles", Gaudi::DataHandle::Reader, this}; - DataHandle m_outputParticles{"SmearedFarForwardParticles", + DataHandle m_outputParticles{"SmearedFarForwardParticles", Gaudi::DataHandle::Writer, this}; - DataHandle m_outputAssocCollection{"MCRecoParticleAssociation", + DataHandle m_outputAssocCollection{"MCRecoParticleAssociation", Gaudi::DataHandle::Writer, this}; Gaudi::Property m_enableZDC{this, "enableZDC", true}; @@ -66,8 +66,8 @@ class SmearedFarForwardParticles : public GaudiAlgorithm { Rndm::Numbers m_gaussDist; - using RecPart = eicd::MutableReconstructedParticle; - using Assoc = eicd::MutableMCRecoParticleAssociation; + using RecPart = edm4eic::MutableReconstructedParticle; + using Assoc = edm4eic::MutableMCRecoParticleAssociation; using RecData = std::pair; public: @@ -160,8 +160,8 @@ class SmearedFarForwardParticles : public GaudiAlgorithm { continue; } // only 0-->4.5 mrad - const double mom_ion_theta = eicd::anglePolar(mom_ion); - const double mom_ion_phi = eicd::angleAzimuthal(mom_ion); + const double mom_ion_theta = edm4eic::anglePolar(mom_ion); + const double mom_ion_phi = edm4eic::angleAzimuthal(mom_ion); if (mom_ion_theta > 4.5 / 1000.) { continue; } @@ -192,7 +192,7 @@ class SmearedFarForwardParticles : public GaudiAlgorithm { const double phis = phi + dphi; const double moms = sqrt(Es * Es - part.getMass() * part.getMass()); // now cast back into float - const auto mom3s_ion = eicd::sphericalToVector(moms, ths, phis); + const auto mom3s_ion = edm4eic::sphericalToVector(moms, ths, phis); const auto mom3s = rotateIonToLabDirection(mom3s_ion); RecPart rec_part; rec_part.setType(kTagZDC); @@ -220,8 +220,8 @@ class SmearedFarForwardParticles : public GaudiAlgorithm { debug() << fmt::format( "Found ZDC particle: {}, Etrue: {}, Emeas: {}, ptrue: {}, pmeas: {}, theta_true: {}, theta_meas: {}", - part.getPDG(), E, rec_part.getEnergy(), part_p_mag, eicd::magnitude(rec_part.getMomentum()), th, - eicd::anglePolar(rec_part.getMomentum())) + part.getPDG(), E, rec_part.getEnergy(), part_p_mag, edm4eic::magnitude(rec_part.getMomentum()), th, + edm4eic::anglePolar(rec_part.getMomentum())) << endmsg; } } @@ -245,7 +245,7 @@ class SmearedFarForwardParticles : public GaudiAlgorithm { } // only 6-->20 mrad const auto mom_ion = removeCrossingAngle(part.getMomentum()); // rotateLabToIonDirection(part.getMomentum()); - const auto mom_ion_theta = eicd::anglePolar(mom_ion); + const auto mom_ion_theta = edm4eic::anglePolar(mom_ion); if (mom_ion_theta < m_thetaMinB0 || mom_ion_theta > m_thetaMaxB0) { continue; } @@ -259,14 +259,14 @@ class SmearedFarForwardParticles : public GaudiAlgorithm { rc.emplace_back(rc_part, assoc); if (msgLevel(MSG::DEBUG)) { const auto& part_p = part.getMomentum(); - const auto part_p_pt = eicd::magnitudeTransverse(part_p); - const auto part_p_mag = eicd::magnitude(part_p); - const auto part_p_theta = eicd::anglePolar(part_p); + const auto part_p_pt = edm4eic::magnitudeTransverse(part_p); + const auto part_p_mag = edm4eic::magnitude(part_p); + const auto part_p_theta = edm4eic::anglePolar(part_p); debug() << fmt::format("Found B0 particle: {}, ptrue: {}, pmeas: {}, pttrue: {}, ptmeas: {}, theta_true: {}, " "theta_meas: {}", - part.getPDG(), part_p_mag, eicd::magnitude(rc_part.momentum()), part_p_pt, - eicd::magnitudeTransverse(rc_part.momentum()), part_p_theta, - eicd::anglePolar(rc_part.momentum())) + part.getPDG(), part_p_mag, edm4eic::magnitude(rc_part.momentum()), part_p_pt, + edm4eic::magnitudeTransverse(rc_part.momentum()), part_p_theta, + edm4eic::anglePolar(rc_part.momentum())) << endmsg; } } @@ -288,7 +288,7 @@ class SmearedFarForwardParticles : public GaudiAlgorithm { continue; } const auto mom_ion = removeCrossingAngle(part.getMomentum()); // rotateLabToIonDirection(part.getMomentum()); - const auto mom_ion_theta = eicd::anglePolar(mom_ion); + const auto mom_ion_theta = edm4eic::anglePolar(mom_ion); if (mom_ion_theta < m_thetaMinRP || mom_ion_theta > m_thetaMaxRP || mom_ion.z < m_pMinRigidityRP * ionBeamEnergy) { continue; @@ -298,14 +298,14 @@ class SmearedFarForwardParticles : public GaudiAlgorithm { rc.emplace_back(rc_part, assoc); if (msgLevel(MSG::DEBUG)) { const auto& part_p = part.getMomentum(); - const auto part_p_pt = eicd::magnitudeTransverse(part_p); - const auto part_p_mag = eicd::magnitude(part_p); - const auto part_p_theta = eicd::anglePolar(part_p); + const auto part_p_pt = edm4eic::magnitudeTransverse(part_p); + const auto part_p_mag = edm4eic::magnitude(part_p); + const auto part_p_theta = edm4eic::anglePolar(part_p); debug() << fmt::format("Found RP particle: {}, ptrue: {}, pmeas: {}, pttrue: {}, ptmeas: {}, theta_true: {}, " "theta_meas: {}", - part.getPDG(), part_p_mag, eicd::magnitude(rc_part.momentum()), part_p_pt, - eicd::magnitudeTransverse(rc_part.momentum()), part_p_theta, - eicd::anglePolar(rc_part.momentum())) + part.getPDG(), part_p_mag, edm4eic::magnitude(rc_part.momentum()), part_p_pt, + edm4eic::magnitudeTransverse(rc_part.momentum()), part_p_theta, + edm4eic::anglePolar(rc_part.momentum())) << endmsg; } } @@ -334,14 +334,14 @@ class SmearedFarForwardParticles : public GaudiAlgorithm { rc.emplace_back(rc_part, assoc); if (msgLevel(MSG::DEBUG)) { const auto& part_p = part.getMomentum(); - const auto part_p_pt = eicd::magnitudeTransverse(part_p); - const auto part_p_mag = eicd::magnitude(part_p); - const auto part_p_theta = eicd::anglePolar(part_p); + const auto part_p_pt = edm4eic::magnitudeTransverse(part_p); + const auto part_p_mag = edm4eic::magnitude(part_p); + const auto part_p_theta = edm4eic::anglePolar(part_p); debug() << fmt::format("Found OMD particle: {}, ptrue: {}, pmeas: {}, pttrue: {}, ptmeas: {}, theta_true: {}, " "theta_meas: {}", - part.getPDG(), part_p_mag, eicd::magnitude(rc_part.momentum()), part_p_pt, - eicd::magnitudeTransverse(rc_part.momentum()), part_p_theta, - eicd::anglePolar(rc_part.momentum())) + part.getPDG(), part_p_mag, edm4eic::magnitude(rc_part.momentum()), part_p_pt, + edm4eic::magnitudeTransverse(rc_part.momentum()), part_p_theta, + edm4eic::anglePolar(rc_part.momentum())) << endmsg; } } @@ -371,7 +371,7 @@ class SmearedFarForwardParticles : public GaudiAlgorithm { // And build our 3-vector const edm4hep::Vector3f psmear_ion{static_cast(pxs), static_cast(pys), static_cast(pzs)}; const auto psmear = rotateIonToLabDirection(psmear_ion); - eicd::MutableReconstructedParticle rec_part; + edm4eic::MutableReconstructedParticle rec_part; rec_part.setType(-1); rec_part.setEnergy(std::hypot(ps, part.getMass())); rec_part.setMomentum({psmear.x, psmear.y, psmear.z}); diff --git a/JugFast/src/components/TruthClustering.cpp b/JugFast/src/components/TruthClustering.cpp index 35ce7ea..7b9e49d 100644 --- a/JugFast/src/components/TruthClustering.cpp +++ b/JugFast/src/components/TruthClustering.cpp @@ -21,10 +21,10 @@ // Event Model related classes #include "edm4hep/MCParticle.h" #include "edm4hep/SimCalorimeterHitCollection.h" -#include "eicd/CalorimeterHitCollection.h" -#include "eicd/ClusterCollection.h" -#include "eicd/ProtoClusterCollection.h" -#include "eicd/RawCalorimeterHitCollection.h" +#include "edm4eic/CalorimeterHitCollection.h" +#include "edm4eic/ClusterCollection.h" +#include "edm4eic/ProtoClusterCollection.h" +#include "edm4eic/RawCalorimeterHitCollection.h" using namespace Gaudi::Units; @@ -36,9 +36,9 @@ namespace Jug::Fast { */ class TruthClustering : public GaudiAlgorithm { private: - DataHandle m_inputHits{"inputHits", Gaudi::DataHandle::Reader, this}; + DataHandle m_inputHits{"inputHits", Gaudi::DataHandle::Reader, this}; DataHandle m_mcHits{"mcHits", Gaudi::DataHandle::Reader, this}; - DataHandle m_outputProtoClusters{"outputProtoClusters", Gaudi::DataHandle::Writer, this}; + DataHandle m_outputProtoClusters{"outputProtoClusters", Gaudi::DataHandle::Writer, this}; public: TruthClustering(const std::string& name, ISvcLocator* svcLoc) diff --git a/JugFast/src/components/TruthEnergyPositionClusterMerger.cpp b/JugFast/src/components/TruthEnergyPositionClusterMerger.cpp index 091cdc2..a743c83 100644 --- a/JugFast/src/components/TruthEnergyPositionClusterMerger.cpp +++ b/JugFast/src/components/TruthEnergyPositionClusterMerger.cpp @@ -16,9 +16,9 @@ // Event Model related classes #include "edm4hep/MCParticleCollection.h" -#include "eicd/ClusterCollection.h" -#include "eicd/MCRecoClusterParticleAssociationCollection.h" -#include +#include "edm4eic/ClusterCollection.h" +#include "edm4eic/MCRecoClusterParticleAssociationCollection.h" +#include using namespace Gaudi::Units; @@ -36,13 +36,13 @@ class TruthEnergyPositionClusterMerger : public GaudiAlgorithm { private: // Input DataHandle m_inputMCParticles{"MCParticles", Gaudi::DataHandle::Reader, this}; - DataHandle m_energyClusters{"EnergyClusters", Gaudi::DataHandle::Reader, this}; - DataHandle m_energyAssociations{"EnergyAssociations", Gaudi::DataHandle::Reader, this}; - DataHandle m_positionClusters{"PositionClusters", Gaudi::DataHandle::Reader, this}; - DataHandle m_positionAssociations{"PositionAssociations", Gaudi::DataHandle::Reader, this}; + DataHandle m_energyClusters{"EnergyClusters", Gaudi::DataHandle::Reader, this}; + DataHandle m_energyAssociations{"EnergyAssociations", Gaudi::DataHandle::Reader, this}; + DataHandle m_positionClusters{"PositionClusters", Gaudi::DataHandle::Reader, this}; + DataHandle m_positionAssociations{"PositionAssociations", Gaudi::DataHandle::Reader, this}; // Output - DataHandle m_outputClusters{"OutputClusters", Gaudi::DataHandle::Writer, this}; - DataHandle m_outputAssociations{"OutputAssociations", Gaudi::DataHandle::Writer, this}; + DataHandle m_outputClusters{"OutputClusters", Gaudi::DataHandle::Writer, this}; + DataHandle m_outputAssociations{"OutputAssociations", Gaudi::DataHandle::Writer, this}; public: TruthEnergyPositionClusterMerger(const std::string& name, ISvcLocator* svcLoc) @@ -121,7 +121,7 @@ class TruthEnergyPositionClusterMerger : public GaudiAlgorithm { } // set association - eicd::MutableMCRecoClusterParticleAssociation clusterassoc; + edm4eic::MutableMCRecoClusterParticleAssociation clusterassoc; clusterassoc.setRecID(new_clus.getObjectID().index); clusterassoc.setSimID(mcID); clusterassoc.setWeight(1.0); @@ -141,7 +141,7 @@ class TruthEnergyPositionClusterMerger : public GaudiAlgorithm { merged_clus.push_back(new_clus); // set association - eicd::MutableMCRecoClusterParticleAssociation clusterassoc; + edm4eic::MutableMCRecoClusterParticleAssociation clusterassoc; clusterassoc.setRecID(new_clus.getObjectID().index); clusterassoc.setSimID(mcID); clusterassoc.setWeight(1.0); @@ -167,7 +167,7 @@ class TruthEnergyPositionClusterMerger : public GaudiAlgorithm { new_clus.setTime(eclus.getTime()); new_clus.setNhits(eclus.getNhits()); // use nominal radius of 110cm, and use start vertex theta and phi - new_clus.setPosition(eicd::sphericalToVector(110.*cm, theta, phi)); + new_clus.setPosition(edm4eic::sphericalToVector(110.*cm, theta, phi)); new_clus.addToClusters(eclus); if (msgLevel(MSG::DEBUG)) { debug() << " --> Processing energy cluster " << eclus.id() << ", mcID: " << mcID << ", energy: " << eclus.getEnergy() @@ -176,7 +176,7 @@ class TruthEnergyPositionClusterMerger : public GaudiAlgorithm { } // set association - eicd::MutableMCRecoClusterParticleAssociation clusterassoc; + edm4eic::MutableMCRecoClusterParticleAssociation clusterassoc; clusterassoc.setRecID(new_clus.getObjectID().index); clusterassoc.setSimID(mcID); clusterassoc.setWeight(1.0); @@ -191,12 +191,12 @@ class TruthEnergyPositionClusterMerger : public GaudiAlgorithm { // get a map of MCParticle index --> cluster // input: cluster_collections --> list of handles to all cluster collections - std::map indexedClusters( - const eicd::ClusterCollection& clusters, - const eicd::MCRecoClusterParticleAssociationCollection& associations + std::map indexedClusters( + const edm4eic::ClusterCollection& clusters, + const edm4eic::MCRecoClusterParticleAssociationCollection& associations ) const { - std::map matched = {}; + std::map matched = {}; for (const auto& cluster : clusters) { int mcID = -1; diff --git a/JugPID/CMakeLists.txt b/JugPID/CMakeLists.txt index 6669b90..c407301 100644 --- a/JugPID/CMakeLists.txt +++ b/JugPID/CMakeLists.txt @@ -14,7 +14,7 @@ gaudi_add_module(JugPIDPlugins JugBase ROOT::Core ROOT::RIO ROOT::Tree EDM4HEP::edm4hep - EICD::eicd + EDM4EIC::edm4eic ) target_include_directories(JugPIDPlugins PUBLIC diff --git a/JugPID/src/components/PhotoRingClusters.cpp b/JugPID/src/components/PhotoRingClusters.cpp index 921c31b..8055405 100644 --- a/JugPID/src/components/PhotoRingClusters.cpp +++ b/JugPID/src/components/PhotoRingClusters.cpp @@ -27,8 +27,8 @@ // Event Model related classes #include "FuzzyKClusters.h" -#include "eicd/PMTHitCollection.h" -#include "eicd/RingImageCollection.h" +#include "edm4eic/PMTHitCollection.h" +#include "edm4eic/RingImageCollection.h" using namespace Gaudi::Units; using namespace Eigen; @@ -41,8 +41,8 @@ namespace Jug::Reco { */ class PhotoRingClusters : public GaudiAlgorithm { private: - DataHandle m_inputHitCollection{"inputHitCollection", Gaudi::DataHandle::Reader, this}; - DataHandle m_outputClusterCollection{"outputClusterCollection", Gaudi::DataHandle::Writer, + DataHandle m_inputHitCollection{"inputHitCollection", Gaudi::DataHandle::Reader, this}; + DataHandle m_outputClusterCollection{"outputClusterCollection", Gaudi::DataHandle::Writer, this}; // @TODO // A more realistic way is to have tracker info as the input to determine how much clusters should be found diff --git a/JugReco/CMakeLists.txt b/JugReco/CMakeLists.txt index d342b85..748d0b2 100644 --- a/JugReco/CMakeLists.txt +++ b/JugReco/CMakeLists.txt @@ -11,10 +11,11 @@ gaudi_add_module(JugRecoPlugins ${JugRecoPlugins_sources} LINK Gaudi::GaudiAlgLib Gaudi::GaudiKernel - JugBase + JugBase JugAlgo + algocore algocalorimetry ROOT::Core ROOT::RIO ROOT::Tree EDM4HEP::edm4hep - EICD::eicd + EDM4EIC::edm4eic DD4hep::DDRec ) diff --git a/JugReco/JugReco/SourceLinks.h b/JugReco/JugReco/SourceLinks.h index 4c5fe8f..85ebb2d 100644 --- a/JugReco/JugReco/SourceLinks.h +++ b/JugReco/JugReco/SourceLinks.h @@ -10,7 +10,7 @@ #include #include "edm4hep/Geant4Particle.h" -#include "eicd/TrackerHitCollection.h" +#include "edm4eic/TrackerHitCollection.h" namespace Jug { @@ -36,7 +36,7 @@ class SourceLink { // need to store pointers to make the object copyable const Acts::Surface* m_surface; //const ActsFatras::Hit* m_truthHit; - const eicd::TrackerHit* m_Hit ; + const edm4eic::TrackerHit* m_Hit ; public: SourceLink(const Acts::Surface& surface, //const ActsFatras::Hit& truthHit, diff --git a/JugReco/src/components/CalorimeterHitReco.cpp b/JugReco/src/components/CalorimeterHitReco.cpp index 7bd08d7..674d7dc 100644 --- a/JugReco/src/components/CalorimeterHitReco.cpp +++ b/JugReco/src/components/CalorimeterHitReco.cpp @@ -26,8 +26,8 @@ #include "JugBase/IGeoSvc.h" // Event Model related classes -#include "eicd/CalorimeterHitCollection.h" -#include "eicd/RawCalorimeterHitCollection.h" +#include "edm4eic/CalorimeterHitCollection.h" +#include "edm4eic/RawCalorimeterHitCollection.h" using namespace Gaudi::Units; @@ -62,9 +62,9 @@ class CalorimeterHitReco : public GaudiAlgorithm { double thresholdADC{0}; double stepTDC{0}; - DataHandle m_inputHitCollection{"inputHitCollection", Gaudi::DataHandle::Reader, + DataHandle m_inputHitCollection{"inputHitCollection", Gaudi::DataHandle::Reader, this}; - DataHandle m_outputHitCollection{"outputHitCollection", Gaudi::DataHandle::Writer, + DataHandle m_outputHitCollection{"outputHitCollection", Gaudi::DataHandle::Writer, this}; // geometry service to get ids, ignored if no names provided @@ -208,13 +208,13 @@ class CalorimeterHitReco : public GaudiAlgorithm { } // create const vectors for passing to hit initializer list - const decltype(eicd::CalorimeterHitData::position) position( + const decltype(edm4eic::CalorimeterHitData::position) position( gpos.x() / m_lUnit, gpos.y() / m_lUnit, gpos.z() / m_lUnit ); - const decltype(eicd::CalorimeterHitData::dimension) dimension( + const decltype(edm4eic::CalorimeterHitData::dimension) dimension( cdim[0] / m_lUnit, cdim[1] / m_lUnit, cdim[2] / m_lUnit ); - const decltype(eicd::CalorimeterHitData::local) local_position( + const decltype(edm4eic::CalorimeterHitData::local) local_position( pos.x() / m_lUnit, pos.y() / m_lUnit, pos.z() / m_lUnit ); diff --git a/JugReco/src/components/CalorimeterHitsEtaPhiProjector.cpp b/JugReco/src/components/CalorimeterHitsEtaPhiProjector.cpp index fe92735..12d414d 100644 --- a/JugReco/src/components/CalorimeterHitsEtaPhiProjector.cpp +++ b/JugReco/src/components/CalorimeterHitsEtaPhiProjector.cpp @@ -32,8 +32,8 @@ #include "JugBase/Utilities/Utils.hpp" // Event Model related classes -#include "eicd/CalorimeterHitCollection.h" -#include "eicd/vector_utils.h" +#include "edm4eic/CalorimeterHitCollection.h" +#include "edm4eic/vector_utils.h" using namespace Gaudi::Units; using Point3D = ROOT::Math::XYZPoint; @@ -60,9 +60,9 @@ namespace Jug::Reco { class CalorimeterHitsEtaPhiProjector : public GaudiAlgorithm { private: Gaudi::Property> u_gridSizes{this, "gridSizes", {0.001, 0.001 * rad}}; - DataHandle m_inputHitCollection{"inputHitCollection", Gaudi::DataHandle::Reader, + DataHandle m_inputHitCollection{"inputHitCollection", Gaudi::DataHandle::Reader, this}; - DataHandle m_outputHitCollection{"outputHitCollection", Gaudi::DataHandle::Writer, + DataHandle m_outputHitCollection{"outputHitCollection", Gaudi::DataHandle::Writer, this}; double gridSizes[2]{0.0, 0.0}; @@ -93,25 +93,25 @@ class CalorimeterHitsEtaPhiProjector : public GaudiAlgorithm { auto& mhits = *m_outputHitCollection.createAndPut(); // container - std::unordered_map, std::vector, pair_hash> merged_hits; + std::unordered_map, std::vector, pair_hash> merged_hits; for (const auto h : *m_inputHitCollection.get()) { auto bins = - std::make_pair(static_cast(pos2bin(eicd::eta(h.getPosition()), gridSizes[0], 0.)), - static_cast(pos2bin(eicd::angleAzimuthal(h.getPosition()), gridSizes[1], 0.))); + std::make_pair(static_cast(pos2bin(edm4eic::eta(h.getPosition()), gridSizes[0], 0.)), + static_cast(pos2bin(edm4eic::angleAzimuthal(h.getPosition()), gridSizes[1], 0.))); merged_hits[bins].push_back(h); } for (const auto& [bins, hits] : merged_hits) { const auto ref = hits.front(); - eicd::MutableCalorimeterHit hit; + edm4eic::MutableCalorimeterHit hit; hit.setCellID(ref.getCellID()); // TODO, we can do timing cut to reject noises hit.setTime(ref.getTime()); - double r = eicd::magnitude(ref.getPosition()); + double r = edm4eic::magnitude(ref.getPosition()); double eta = bin2pos(bins.first, gridSizes[0], 0.); double phi = bin2pos(bins.second, gridSizes[1], 1.); - hit.setPosition(eicd::sphericalToVector(r, eicd::etaToAngle(eta), phi)); + hit.setPosition(edm4eic::sphericalToVector(r, edm4eic::etaToAngle(eta), phi)); hit.setDimension({static_cast(gridSizes[0]), static_cast(gridSizes[1]), 0.}); // merge energy hit.setEnergy(0.); diff --git a/JugReco/src/components/CalorimeterHitsMerger.cpp b/JugReco/src/components/CalorimeterHitsMerger.cpp index 3b7c0b6..db0becc 100644 --- a/JugReco/src/components/CalorimeterHitsMerger.cpp +++ b/JugReco/src/components/CalorimeterHitsMerger.cpp @@ -32,7 +32,7 @@ #include "JugBase/IGeoSvc.h" // Event Model related classes -#include "eicd/CalorimeterHitCollection.h" +#include "edm4eic/CalorimeterHitCollection.h" using namespace Gaudi::Units; @@ -53,8 +53,8 @@ class CalorimeterHitsMerger : public GaudiAlgorithm { Gaudi::Property> u_fields{this, "fields", {"layer"}}; // reference field numbers to locate position for each merged hits group Gaudi::Property> u_refs{this, "fieldRefNumbers", {}}; - DataHandle m_inputHitCollection{"inputHitCollection", Gaudi::DataHandle::Reader, this}; - DataHandle m_outputHitCollection{"outputHitCollection", Gaudi::DataHandle::Writer, + DataHandle m_inputHitCollection{"inputHitCollection", Gaudi::DataHandle::Reader, this}; + DataHandle m_outputHitCollection{"outputHitCollection", Gaudi::DataHandle::Writer, this}; SmartIF m_geoSvc; @@ -111,7 +111,7 @@ class CalorimeterHitsMerger : public GaudiAlgorithm { auto& outputs = *m_outputHitCollection.createAndPut(); // find the hits that belong to the same group (for merging) - std::unordered_map> merge_map; + std::unordered_map> merge_map; for (const auto& h : inputs) { int64_t id = h.getCellID() & id_mask; // use the reference field position @@ -162,15 +162,15 @@ class CalorimeterHitsMerger : public GaudiAlgorithm { const auto& href = hits.front(); // create const vectors for passing to hit initializer list - const decltype(eicd::CalorimeterHitData::position) position( + const decltype(edm4eic::CalorimeterHitData::position) position( gpos.x() / dd4hep::mm, gpos.y() / dd4hep::mm, gpos.z() / dd4hep::mm ); - const decltype(eicd::CalorimeterHitData::local) local( + const decltype(edm4eic::CalorimeterHitData::local) local( pos.x(), pos.y(), pos.z() ); outputs.push_back( - eicd::CalorimeterHit{href.getCellID(), + edm4eic::CalorimeterHit{href.getCellID(), energy, energyError, time, diff --git a/JugReco/src/components/CalorimeterIslandCluster.cpp b/JugReco/src/components/CalorimeterIslandCluster.cpp index 194f647..3a3ab6d 100644 --- a/JugReco/src/components/CalorimeterIslandCluster.cpp +++ b/JugReco/src/components/CalorimeterIslandCluster.cpp @@ -33,64 +33,84 @@ #include "JugBase/IGeoSvc.h" // Event Model related classes -#include "eicd/CalorimeterHitCollection.h" -#include "eicd/ClusterCollection.h" -#include "eicd/ProtoClusterCollection.h" -#include "eicd/Vector2f.h" -#include "eicd/Vector3f.h" -#include "eicd/vector_utils.h" +#include "edm4eic/CalorimeterHitCollection.h" +#include "edm4eic/ClusterCollection.h" +#include "edm4eic/ProtoClusterCollection.h" +#include "edm4eic/vector_utils.h" +#include "edm4hep/Vector3f.h" +#include "edm4hep/Vector2f.h" + +#if defined __has_include +# if __has_include ("edm4eic/Vector3f.h") +# include "edm4eic/Vector3f.h" +# endif +# if __has_include ("edm4eic/Vector2f.h") +# include "edm4eic/Vector2f.h" +# endif +#endif + +namespace edm4eic { + class Vector2f; + class Vector3f; +} using namespace Gaudi::Units; namespace { -using CaloHit = eicd::CalorimeterHit; -using CaloHitCollection = eicd::CalorimeterHitCollection; +using CaloHit = edm4eic::CalorimeterHit; +using CaloHitCollection = edm4eic::CalorimeterHitCollection; + +using Vector2f = std::conditional_t< + std::is_same_v, + edm4hep::Vector2f, + edm4eic::Vector2f +>; // helper functions to get distance between hits -static eicd::Vector2f localDistXY(const CaloHit& h1, const CaloHit& h2) { +static Vector2f localDistXY(const CaloHit& h1, const CaloHit& h2) { const auto delta = h1.getLocal() - h2.getLocal(); return {delta.x, delta.y}; } -static eicd::Vector2f localDistXZ(const CaloHit& h1, const CaloHit& h2) { +static Vector2f localDistXZ(const CaloHit& h1, const CaloHit& h2) { const auto delta = h1.getLocal() - h2.getLocal(); return {delta.x, delta.z}; } -static eicd::Vector2f localDistYZ(const CaloHit& h1, const CaloHit& h2) { +static Vector2f localDistYZ(const CaloHit& h1, const CaloHit& h2) { const auto delta = h1.getLocal() - h2.getLocal(); return {delta.y, delta.z}; } -static eicd::Vector2f dimScaledLocalDistXY(const CaloHit& h1, const CaloHit& h2) { +static Vector2f dimScaledLocalDistXY(const CaloHit& h1, const CaloHit& h2) { const auto delta = h1.getLocal() - h2.getLocal(); const auto dimsum = h1.getDimension() + h2.getDimension(); return {2 * delta.x / dimsum.x, 2 * delta.y / dimsum.y}; } -static eicd::Vector2f globalDistRPhi(const CaloHit& h1, const CaloHit& h2) { - using vector_type = decltype(eicd::Vector2f::a); +static Vector2f globalDistRPhi(const CaloHit& h1, const CaloHit& h2) { + using vector_type = decltype(Vector2f::a); return { static_cast( - eicd::magnitude(h1.getPosition()) - eicd::magnitude(h2.getPosition()) + edm4eic::magnitude(h1.getPosition()) - edm4eic::magnitude(h2.getPosition()) ), static_cast( - eicd::angleAzimuthal(h1.getPosition()) - eicd::angleAzimuthal(h2.getPosition()) + edm4eic::angleAzimuthal(h1.getPosition()) - edm4eic::angleAzimuthal(h2.getPosition()) ) }; } -static eicd::Vector2f globalDistEtaPhi(const CaloHit& h1, +static Vector2f globalDistEtaPhi(const CaloHit& h1, const CaloHit& h2) { - using vector_type = decltype(eicd::Vector2f::a); + using vector_type = decltype(Vector2f::a); return { static_cast( - eicd::eta(h1.getPosition()) - eicd::eta(h2.getPosition()) + edm4eic::eta(h1.getPosition()) - edm4eic::eta(h2.getPosition()) ), static_cast( - eicd::angleAzimuthal(h1.getPosition()) - eicd::angleAzimuthal(h2.getPosition()) + edm4eic::angleAzimuthal(h1.getPosition()) - edm4eic::angleAzimuthal(h2.getPosition()) ) }; } // name: {method, units} static std::map, std::vector>> + std::tuple, std::vector>> distMethods{ {"localDistXY", {localDistXY, {mm, mm}}}, {"localDistXZ", {localDistXZ, {mm, mm}}}, {"localDistYZ", {localDistYZ, {mm, mm}}}, {"dimScaledLocalDistXY", {dimScaledLocalDistXY, {1., 1.}}}, @@ -118,7 +138,7 @@ class CalorimeterIslandCluster : public GaudiAlgorithm { Gaudi::Property m_minClusterHitEdep{this, "minClusterHitEdep", 0.}; Gaudi::Property m_minClusterCenterEdep{this, "minClusterCenterEdep", 50.0 * MeV}; DataHandle m_inputHitCollection{"inputHitCollection", Gaudi::DataHandle::Reader, this}; - DataHandle m_outputProtoCollection{"outputProtoClusterCollection", + DataHandle m_outputProtoCollection{"outputProtoClusterCollection", Gaudi::DataHandle::Writer, this}; // neighbour checking distances @@ -130,7 +150,7 @@ class CalorimeterIslandCluster : public GaudiAlgorithm { Gaudi::Property> u_globalDistEtaPhi{this, "globalDistEtaPhi", {}}; Gaudi::Property> u_dimScaledLocalDistXY{this, "dimScaledLocalDistXY", {1.8, 1.8}}; // neighbor checking function - std::function hitsDist; + std::function hitsDist; // unitless counterparts of the input parameters double minClusterHitEdep{0}, minClusterCenterEdep{0}, sectorDist{0}; @@ -253,7 +273,7 @@ class CalorimeterIslandCluster : public GaudiAlgorithm { // different sector, local coordinates do not work, using global coordinates } else { // sector may have rotation (barrel), so z is included - return (eicd::magnitude(h1.getPosition() - h2.getPosition()) <= sectorDist); + return (edm4eic::magnitude(h1.getPosition() - h2.getPosition()) <= sectorDist); } } @@ -337,7 +357,7 @@ class CalorimeterIslandCluster : public GaudiAlgorithm { // split a group of hits according to the local maxima void split_group(std::vector>& group, const std::vector& maxima, - eicd::ProtoClusterCollection& proto) const { + edm4eic::ProtoClusterCollection& proto) const { // special cases if (maxima.empty()) { if (msgLevel(MSG::VERBOSE)) { @@ -345,7 +365,7 @@ class CalorimeterIslandCluster : public GaudiAlgorithm { } return; } else if (maxima.size() == 1) { - eicd::MutableProtoCluster pcl; + edm4eic::MutableProtoCluster pcl; for (auto& [idx, hit] : group) { pcl.addToHits(hit); pcl.addToWeights(1.); @@ -360,7 +380,7 @@ class CalorimeterIslandCluster : public GaudiAlgorithm { // split between maxima // TODO, here we can implement iterations with profile, or even ML for better splits std::vector weights(maxima.size(), 1.); - std::vector pcls; + std::vector pcls; for (size_t k = 0; k < maxima.size(); ++k) { pcls.emplace_back(); } @@ -372,7 +392,7 @@ class CalorimeterIslandCluster : public GaudiAlgorithm { for (const auto& chit : maxima) { double dist_ref = chit.getDimension().x; double energy = chit.getEnergy(); - double dist = eicd::magnitude(hitsDist(chit, hit)); + double dist = edm4eic::magnitude(hitsDist(chit, hit)); weights[j] = std::exp(-dist / dist_ref) * energy; j += 1; } diff --git a/JugReco/src/components/ClusterRecoCoG.cpp b/JugReco/src/components/ClusterRecoCoG.cpp index 8b60f02..2839ca6 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 "eicd/ClusterCollection.h" -#include "eicd/MCRecoClusterParticleAssociationCollection.h" -#include "eicd/ProtoClusterCollection.h" -#include "eicd/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 - eicd::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() << ", " << eicd::magnitude(mcp.getMomentum()) << endmsg; - } - - // set association - eicd::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: - eicd::MutableCluster reconstruct(const eicd::ProtoCluster& pcl) const { - eicd::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 = eicd::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 = (eicd::eta(cl.getPosition()) > maxHitEta); - const bool underflow = (eicd::eta(cl.getPosition()) < minHitEta); - if (overflow || underflow) { - const double newEta = overflow ? maxHitEta : minHitEta; - const double newTheta = eicd::etaToAngle(newEta); - const double newR = eicd::magnitude(cl.getPosition()); - const double newPhi = eicd::angleAzimuthal(cl.getPosition()); - cl.setPosition(eicd::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(eicd::anglePolar(cl.getPosition())); - cl.setIntrinsicPhi(eicd::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/JugReco/src/components/EnergyPositionClusterMerger.cpp b/JugReco/src/components/EnergyPositionClusterMerger.cpp index ff3aa38..957611e 100644 --- a/JugReco/src/components/EnergyPositionClusterMerger.cpp +++ b/JugReco/src/components/EnergyPositionClusterMerger.cpp @@ -15,8 +15,8 @@ #include "JugBase/DataHandle.h" // Event Model related classes -#include "eicd/ClusterCollection.h" -#include "eicd/vector_utils.h" +#include "edm4eic/ClusterCollection.h" +#include "edm4eic/vector_utils.h" using namespace Gaudi::Units; @@ -37,10 +37,10 @@ namespace Jug::Reco { class EnergyPositionClusterMerger : public GaudiAlgorithm { private: // Input - DataHandle m_energyClusters{"energyClusters", Gaudi::DataHandle::Reader, this}; - DataHandle m_positionClusters{"positionClusters", Gaudi::DataHandle::Reader, this}; + DataHandle m_energyClusters{"energyClusters", Gaudi::DataHandle::Reader, this}; + DataHandle m_positionClusters{"positionClusters", Gaudi::DataHandle::Reader, this}; // Output - DataHandle m_outputClusters{"outputClusters", Gaudi::DataHandle::Writer, this}; + DataHandle m_outputClusters{"outputClusters", Gaudi::DataHandle::Writer, this}; // Negative values mean the tolerance check is disabled Gaudi::Property m_zToleranceUnits{this, "zTolerance", -1 * cm}; Gaudi::Property m_phiToleranceUnits{this, "phiTolerance", 20 * degree}; @@ -83,7 +83,7 @@ class EnergyPositionClusterMerger : public GaudiAlgorithm { const auto& ec = e_clus[ie]; // 1. stop if not within tolerance // (make sure to handle rollover of phi properly) - double dphi = eicd::angleAzimuthal(pc.getPosition()) - eicd::angleAzimuthal(ec.getPosition()); + double dphi = edm4eic::angleAzimuthal(pc.getPosition()) - edm4eic::angleAzimuthal(ec.getPosition()); if (std::abs(dphi) > M_PI) { dphi = std::abs(dphi) - M_PI; } @@ -120,13 +120,13 @@ class EnergyPositionClusterMerger : public GaudiAlgorithm { if (msgLevel(MSG::DEBUG)) { debug() << fmt::format("Matched position cluster {} with energy cluster {}\n", pc.id(), ec.id()) << endmsg; debug() << fmt::format(" - Position cluster: (E: {}, phi: {}, z: {})", pc.getEnergy(), - eicd::angleAzimuthal(pc.getPosition()), pc.getPosition().z) + edm4eic::angleAzimuthal(pc.getPosition()), pc.getPosition().z) << endmsg; debug() << fmt::format(" - Energy cluster: (E: {}, phi: {}, z: {})", ec.getEnergy(), - eicd::angleAzimuthal(ec.getPosition()), ec.getPosition().z) + edm4eic::angleAzimuthal(ec.getPosition()), ec.getPosition().z) << endmsg; debug() << fmt::format(" ---> Merged cluster: (E: {}, phi: {}, z: {})", new_clus.getEnergy(), - eicd::angleAzimuthal(new_clus.getPosition()), new_clus.getPosition().z) + edm4eic::angleAzimuthal(new_clus.getPosition()), new_clus.getPosition().z) << endmsg; } } diff --git a/JugReco/src/components/FarForwardParticles.cpp b/JugReco/src/components/FarForwardParticles.cpp index 828ee36..9922644 100644 --- a/JugReco/src/components/FarForwardParticles.cpp +++ b/JugReco/src/components/FarForwardParticles.cpp @@ -19,16 +19,16 @@ #include "JugBase/IGeoSvc.h" // Event Model related classes -#include "eicd/ReconstructedParticleCollection.h" -#include "eicd/TrackerHitCollection.h" -#include +#include "edm4eic/ReconstructedParticleCollection.h" +#include "edm4eic/TrackerHitCollection.h" +#include namespace Jug::Reco { class FarForwardParticles : public GaudiAlgorithm { private: - DataHandle m_inputHitCollection{"FarForwardTrackerHits", Gaudi::DataHandle::Reader, this}; - DataHandle m_outputParticles{"outputParticles", Gaudi::DataHandle::Writer, + DataHandle m_inputHitCollection{"FarForwardTrackerHits", Gaudi::DataHandle::Reader, this}; + DataHandle m_outputParticles{"outputParticles", Gaudi::DataHandle::Writer, this}; //----- Define constants here ------ @@ -152,7 +152,7 @@ class FarForwardParticles : public GaudiAlgorithm { } StatusCode execute() override { - const eicd::TrackerHitCollection* rawhits = m_inputHitCollection.get(); + const edm4eic::TrackerHitCollection* rawhits = m_inputHitCollection.get(); auto& rc = *(m_outputParticles.createAndPut()); auto converter = m_geoSvc->cellIDPositionConverter(); @@ -255,10 +255,10 @@ class FarForwardParticles : public GaudiAlgorithm { //----- end RP reconstruction code ------ - eicd::MutableReconstructedParticle rpTrack; + edm4eic::MutableReconstructedParticle rpTrack; rpTrack.setType(0); rpTrack.setMomentum({prec}); - rpTrack.setEnergy(std::hypot(eicd::magnitude(rpTrack.getMomentum()), .938272)); + rpTrack.setEnergy(std::hypot(edm4eic::magnitude(rpTrack.getMomentum()), .938272)); rpTrack.setReferencePoint({0, 0, 0}); rpTrack.setCharge(1); rpTrack.setMass(.938272); diff --git a/JugReco/src/components/FarForwardParticlesOMD.cpp b/JugReco/src/components/FarForwardParticlesOMD.cpp index 3a14f8c..49e5514 100644 --- a/JugReco/src/components/FarForwardParticlesOMD.cpp +++ b/JugReco/src/components/FarForwardParticlesOMD.cpp @@ -14,16 +14,16 @@ #include "JugBase/DataHandle.h" // Event Model related classes -#include "eicd/ReconstructedParticleCollection.h" -#include "eicd/TrackerHitCollection.h" -#include +#include "edm4eic/ReconstructedParticleCollection.h" +#include "edm4eic/TrackerHitCollection.h" +#include namespace Jug::Reco { class FarForwardParticlesOMD : public GaudiAlgorithm { private: - DataHandle m_inputHitCollection{"FarForwardTrackerHits", Gaudi::DataHandle::Reader, this}; - DataHandle m_outputParticles{"outputParticles", Gaudi::DataHandle::Writer, + DataHandle m_inputHitCollection{"FarForwardTrackerHits", Gaudi::DataHandle::Reader, this}; + DataHandle m_outputParticles{"outputParticles", Gaudi::DataHandle::Writer, this}; //----- Define constants here ------ @@ -75,7 +75,7 @@ class FarForwardParticlesOMD : public GaudiAlgorithm { } StatusCode execute() override { - const eicd::TrackerHitCollection* rawhits = m_inputHitCollection.get(); + const edm4eic::TrackerHitCollection* rawhits = m_inputHitCollection.get(); auto& rc = *(m_outputParticles.createAndPut()); // for (const auto& part : mc) { @@ -164,10 +164,10 @@ class FarForwardParticlesOMD : public GaudiAlgorithm { //----- end RP reconstruction code ------ - eicd::MutableReconstructedParticle rpTrack; + edm4eic::MutableReconstructedParticle rpTrack; rpTrack.setType(0); rpTrack.setMomentum({prec}); - rpTrack.setEnergy(std::hypot(eicd::magnitude(rpTrack.getMomentum()), .938272)); + rpTrack.setEnergy(std::hypot(edm4eic::magnitude(rpTrack.getMomentum()), .938272)); rpTrack.setReferencePoint({0, 0, 0}); rpTrack.setCharge(1); rpTrack.setMass(.938272); diff --git a/JugReco/src/components/ImagingClusterReco.cpp b/JugReco/src/components/ImagingClusterReco.cpp index 5d78e7d..0e85ec9 100644 --- a/JugReco/src/components/ImagingClusterReco.cpp +++ b/JugReco/src/components/ImagingClusterReco.cpp @@ -31,11 +31,11 @@ // Event Model related classes #include "edm4hep/MCParticleCollection.h" #include "edm4hep/SimCalorimeterHitCollection.h" -#include "eicd/CalorimeterHitCollection.h" -#include "eicd/ClusterCollection.h" -#include "eicd/MCRecoClusterParticleAssociationCollection.h" -#include "eicd/ProtoClusterCollection.h" -#include "eicd/vector_utils.h" +#include "edm4eic/CalorimeterHitCollection.h" +#include "edm4eic/ClusterCollection.h" +#include "edm4eic/MCRecoClusterParticleAssociationCollection.h" +#include "edm4eic/ProtoClusterCollection.h" +#include "edm4eic/vector_utils.h" using namespace Gaudi::Units; using namespace Eigen; @@ -53,9 +53,9 @@ class ImagingClusterReco : public GaudiAlgorithm { private: Gaudi::Property m_trackStopLayer{this, "trackStopLayer", 9}; - DataHandle m_inputProtoClusters{"inputProtoClusters", Gaudi::DataHandle::Reader, this}; - DataHandle m_outputLayers{"outputLayers", Gaudi::DataHandle::Writer, this}; - DataHandle m_outputClusters{"outputClusters", Gaudi::DataHandle::Reader, this}; + DataHandle m_inputProtoClusters{"inputProtoClusters", Gaudi::DataHandle::Reader, this}; + DataHandle m_outputLayers{"outputLayers", Gaudi::DataHandle::Writer, this}; + DataHandle m_outputClusters{"outputClusters", Gaudi::DataHandle::Reader, this}; // Collection for MC hits when running on MC Gaudi::Property m_mcHits{this, "mcHits", ""}; @@ -65,7 +65,7 @@ class ImagingClusterReco : public GaudiAlgorithm { // Collection for associations when running on MC Gaudi::Property m_outputAssociations{this, "outputAssociations", ""}; // Optional handle to MC hits - std::unique_ptr> m_outputAssociations_ptr; + std::unique_ptr> m_outputAssociations_ptr; public: ImagingClusterReco(const std::string& name, ISvcLocator* svcLoc) : GaudiAlgorithm(name, svcLoc) { @@ -89,7 +89,7 @@ class ImagingClusterReco : public GaudiAlgorithm { // Initialize the optional association collection if requested if (m_outputAssociations != "") { m_outputAssociations_ptr = - std::make_unique>(m_outputAssociations, Gaudi::DataHandle::Writer, + std::make_unique>(m_outputAssociations, Gaudi::DataHandle::Writer, this); } @@ -110,7 +110,7 @@ class ImagingClusterReco : public GaudiAlgorithm { } // Optional output associations - eicd::MCRecoClusterParticleAssociationCollection* associations = nullptr; + edm4eic::MCRecoClusterParticleAssociationCollection* associations = nullptr; if (m_outputAssociations_ptr) { associations = m_outputAssociations_ptr->createAndPut(); } @@ -170,7 +170,7 @@ class ImagingClusterReco : public GaudiAlgorithm { const auto& mcp = mchit->getContributions(0).getParticle(); // set association - eicd::MutableMCRecoClusterParticleAssociation clusterassoc; + edm4eic::MutableMCRecoClusterParticleAssociation clusterassoc; clusterassoc.setRecID(cl.getObjectID().index); clusterassoc.setSimID(mcp.getObjectID().index); clusterassoc.setWeight(1.0); @@ -197,11 +197,11 @@ class ImagingClusterReco : public GaudiAlgorithm { private: template static inline T pow2(const T& x) { return x * x; } - static std::vector reconstruct_cluster_layers(const eicd::ProtoCluster& pcl) { + static std::vector reconstruct_cluster_layers(const edm4eic::ProtoCluster& pcl) { const auto& hits = pcl.getHits(); const auto& weights = pcl.getWeights(); // using map to have hits sorted by layer - std::map>> layer_map; + std::map>> layer_map; for (unsigned i = 0; i < hits.size(); ++i) { const auto hit = hits[i]; auto lid = hit.getLayer(); @@ -212,7 +212,7 @@ class ImagingClusterReco : public GaudiAlgorithm { } // create layers - std::vector cl_layers; + std::vector cl_layers; for (const auto& [lid, layer_hits] : layer_map) { auto layer = reconstruct_layer(layer_hits); cl_layers.push_back(layer); @@ -220,8 +220,8 @@ class ImagingClusterReco : public GaudiAlgorithm { return cl_layers; } - static eicd::Cluster reconstruct_layer(const std::vector>& hits) { - eicd::MutableCluster layer; + static edm4eic::Cluster reconstruct_layer(const std::vector>& hits) { + edm4eic::MutableCluster layer; layer.setType(ClusterType::kClusterSlice); // Calculate averages double energy{0}; @@ -251,7 +251,7 @@ class ImagingClusterReco : public GaudiAlgorithm { // Calculate radius as the standard deviation of the hits versus the cluster center double radius = 0.; for (const auto& [hit, weight] : hits) { - radius += std::pow(eicd::magnitude(hit.getPosition() - layer.getPosition()), 2); + radius += std::pow(edm4eic::magnitude(hit.getPosition() - layer.getPosition()), 2); } layer.addToShapeParameters(std::sqrt(radius / layer.getNhits())); // TODO Skewedness @@ -259,8 +259,8 @@ class ImagingClusterReco : public GaudiAlgorithm { return layer; } - eicd::MutableCluster reconstruct_cluster(const eicd::ProtoCluster& pcl) { - eicd::MutableCluster cluster; + edm4eic::MutableCluster reconstruct_cluster(const edm4eic::ProtoCluster& pcl) { + edm4eic::MutableCluster cluster; const auto& hits = pcl.getHits(); const auto& weights = pcl.getWeights(); @@ -282,9 +282,9 @@ class ImagingClusterReco : public GaudiAlgorithm { const double energyWeight = hit.getEnergy() * weight; time += hit.getTime() * energyWeight; timeError += std::pow(hit.getTimeError() * energyWeight, 2); - meta += eicd::eta(hit.getPosition()) * energyWeight; - mphi += eicd::angleAzimuthal(hit.getPosition()) * energyWeight; - r = std::min(eicd::magnitude(hit.getPosition()), r); + meta += edm4eic::eta(hit.getPosition()) * energyWeight; + mphi += edm4eic::angleAzimuthal(hit.getPosition()) * energyWeight; + r = std::min(edm4eic::magnitude(hit.getPosition()), r); cluster.addToHits(hit); } cluster.setEnergy(energy); @@ -292,13 +292,13 @@ class ImagingClusterReco : public GaudiAlgorithm { cluster.setTime(time / energy); cluster.setTimeError(std::sqrt(timeError) / energy); cluster.setNhits(hits.size()); - cluster.setPosition(eicd::sphericalToVector(r, eicd::etaToAngle(meta / energy), mphi / energy)); + cluster.setPosition(edm4eic::sphericalToVector(r, edm4eic::etaToAngle(meta / energy), mphi / energy)); // shower radius estimate (eta-phi plane) double radius = 0.; for (const auto& hit : hits) { - radius += pow2(eicd::eta(hit.getPosition()) - eicd::eta(cluster.getPosition())) + - pow2(eicd::angleAzimuthal(hit.getPosition()) - eicd::angleAzimuthal(cluster.getPosition())); + radius += pow2(edm4eic::eta(hit.getPosition()) - edm4eic::eta(cluster.getPosition())) + + pow2(edm4eic::angleAzimuthal(hit.getPosition()) - edm4eic::angleAzimuthal(cluster.getPosition())); } cluster.addToShapeParameters(std::sqrt(radius / cluster.getNhits())); // Skewedness not calculated TODO @@ -313,9 +313,9 @@ class ImagingClusterReco : public GaudiAlgorithm { return cluster; } - std::pair fit_track(const std::vector& layers) const { + std::pair fit_track(const std::vector& layers) const { int nrows = 0; - eicd::Vector3f mean_pos{0, 0, 0}; + decltype(edm4eic::ClusterData::position) mean_pos{0, 0, 0}; for (const auto& layer : layers) { if ((layer.getNhits() > 0) && (layer.getHits(0).getLayer() <= m_trackStopLayer)) { mean_pos = mean_pos + layer.getPosition(); diff --git a/JugReco/src/components/ImagingPixelDataCombiner.cpp b/JugReco/src/components/ImagingPixelDataCombiner.cpp index dc1c139..4b71e49 100644 --- a/JugReco/src/components/ImagingPixelDataCombiner.cpp +++ b/JugReco/src/components/ImagingPixelDataCombiner.cpp @@ -28,8 +28,8 @@ #include "JugBase/Utilities/Utils.hpp" // Event Model related classes -#include "eicd/CalorimeterHitCollection.h" -#include "eicd/vector_utils.h" +#include "edm4eic/CalorimeterHitCollection.h" +#include "edm4eic/vector_utils.h" using namespace Gaudi::Units; @@ -48,9 +48,9 @@ class ImagingPixelDataCombiner : public GaudiAlgorithm { private: Gaudi::Property m_layerIncrement{this, "layerIncrement", 0}; Gaudi::Property m_rule{this, "rule", "concatenate"}; - DataHandle m_inputHits1{"inputHits1", Gaudi::DataHandle::Reader, this}; - DataHandle m_inputHits2{"inputHits2", Gaudi::DataHandle::Reader, this}; - DataHandle m_outputHits{"outputHits", Gaudi::DataHandle::Writer, this}; + DataHandle m_inputHits1{"inputHits1", Gaudi::DataHandle::Reader, this}; + DataHandle m_inputHits2{"inputHits2", Gaudi::DataHandle::Reader, this}; + DataHandle m_outputHits{"outputHits", Gaudi::DataHandle::Writer, this}; std::vector supported_rules{"concatenate", "interlayer"}; public: @@ -80,7 +80,7 @@ class ImagingPixelDataCombiner : public GaudiAlgorithm { // input collections const auto* const hits1 = m_inputHits1.get(); const auto* const hits2 = m_inputHits2.get(); - std::vector inputs{hits1, hits2}; + std::vector inputs{hits1, hits2}; // Create output collections auto* mhits = m_outputHits.createAndPut(); @@ -89,7 +89,7 @@ class ImagingPixelDataCombiner : public GaudiAlgorithm { for (int i = 0; i < (int)inputs.size(); ++i) { const auto* const coll = inputs[i]; for (auto hit : *coll) { - eicd::CalorimeterHit h2{ + edm4eic::CalorimeterHit h2{ hit.getCellID(), hit.getEnergy(), hit.getEnergyError(), hit.getTime(), hit.getTimeError(), hit.getPosition(), hit.getDimension(), hit.getLayer() + m_layerIncrement * i, hit.getSector(), hit.getLocal(), @@ -140,7 +140,7 @@ class ImagingPixelDataCombiner : public GaudiAlgorithm { } // push hit, increment of index - eicd::CalorimeterHit h2{ + edm4eic::CalorimeterHit h2{ hit.getCellID(), hit.getEnergy(), hit.getEnergyError(), hit.getTime(), hit.getTimeError(), hit.getPosition(), hit.getDimension(), hit.getLayer() + m_layerIncrement * curr_coll, hit.getSector(), hit.getLocal()}; diff --git a/JugReco/src/components/ImagingPixelDataSorter.cpp b/JugReco/src/components/ImagingPixelDataSorter.cpp index 84cea6c..83484b7 100644 --- a/JugReco/src/components/ImagingPixelDataSorter.cpp +++ b/JugReco/src/components/ImagingPixelDataSorter.cpp @@ -27,9 +27,9 @@ #include "JugBase/Utilities/Utils.hpp" // Event Model related classes -#include -#include "eicd/CalorimeterHit.h" -#include "eicd/CalorimeterHitCollection.h" +#include +#include "edm4eic/CalorimeterHit.h" +#include "edm4eic/CalorimeterHitCollection.h" using namespace Gaudi::Units; @@ -48,9 +48,9 @@ namespace Jug::Reco { private: Gaudi::Property m_nLayers{this, "numberOfLayers", 9}; Gaudi::Property m_nHits{this, "numberOfHits", 50}; - DataHandle m_inputHitCollection{"inputHitCollection", + DataHandle m_inputHitCollection{"inputHitCollection", Gaudi::DataHandle::Reader, this}; - DataHandle m_outputHitCollection{"outputHitCollection", + DataHandle m_outputHitCollection{"outputHitCollection", Gaudi::DataHandle::Writer, this}; public: @@ -78,7 +78,7 @@ namespace Jug::Reco { auto& mhits = *m_outputHitCollection.createAndPut(); // group the hits by layer - std::vector> layer_hits; + std::vector> layer_hits; layer_hits.resize(m_nLayers); for (const auto& h : hits) { auto k = h.getLayer(); @@ -90,7 +90,7 @@ namespace Jug::Reco { // sort by energy for (auto &layer : layer_hits) { std::sort(layer.begin(), layer.end(), - [] (const eicd::CalorimeterHit &h1, const eicd::CalorimeterHit &h2) { + [] (const edm4eic::CalorimeterHit &h1, const edm4eic::CalorimeterHit &h2) { return h1.getEnergy() > h2.getEnergy(); }); } diff --git a/JugReco/src/components/ImagingPixelMerger.cpp b/JugReco/src/components/ImagingPixelMerger.cpp index f641a94..c81219a 100644 --- a/JugReco/src/components/ImagingPixelMerger.cpp +++ b/JugReco/src/components/ImagingPixelMerger.cpp @@ -31,8 +31,8 @@ #include "JugBase/Utilities/Utils.hpp" // Event Model related classes -#include "eicd/CalorimeterHitCollection.h" -#include +#include "edm4eic/CalorimeterHitCollection.h" +#include using namespace Gaudi::Units; @@ -56,8 +56,8 @@ class ImagingPixelMerger : public GaudiAlgorithm { private: Gaudi::Property m_etaSize{this, "etaSize", 0.001}; Gaudi::Property m_phiSize{this, "phiSize", 0.001}; - DataHandle m_inputHits{"inputHits", Gaudi::DataHandle::Reader, this}; - DataHandle m_outputHits{"outputHits", Gaudi::DataHandle::Writer, this}; + DataHandle m_inputHits{"inputHits", Gaudi::DataHandle::Reader, this}; + DataHandle m_outputHits{"outputHits", Gaudi::DataHandle::Writer, this}; public: ImagingPixelMerger(const std::string& name, ISvcLocator* svcLoc) : GaudiAlgorithm(name, svcLoc) { @@ -102,9 +102,9 @@ class ImagingPixelMerger : public GaudiAlgorithm { const auto& pos = h.getPosition(); // cylindrical r - const float rc = eicd::magnitudeTransverse(pos); - const double eta = eicd::eta(pos); - const double phi = eicd::angleAzimuthal(pos); + const float rc = edm4eic::magnitudeTransverse(pos); + const double eta = edm4eic::eta(pos); + const double phi = edm4eic::angleAzimuthal(pos); const auto grid = std::pair{pos2grid(eta, m_etaSize), pos2grid(phi, m_phiSize)}; auto it = layer.find(grid); @@ -132,10 +132,10 @@ class ImagingPixelMerger : public GaudiAlgorithm { for (const auto& [grid, data] : layer) { const double eta = grid2pos(grid.first, m_etaSize); const double phi = grid2pos(grid.second, m_phiSize); - const double theta = eicd::etaToAngle(eta); + const double theta = edm4eic::etaToAngle(eta); const double z = cotan(theta) * data.rc; const float r = std::hypot(data.rc, z); - const auto pos = eicd::sphericalToVector(r, theta, phi); + const auto pos = edm4eic::sphericalToVector(r, theta, phi); auto oh = ohits.create(); oh.setEnergy(data.energy); oh.setEnergyError(std::sqrt(data.energyError)); diff --git a/JugReco/src/components/ImagingPixelReco.cpp b/JugReco/src/components/ImagingPixelReco.cpp index 5ab8318..f815f07 100644 --- a/JugReco/src/components/ImagingPixelReco.cpp +++ b/JugReco/src/components/ImagingPixelReco.cpp @@ -24,8 +24,8 @@ #include "JugBase/IGeoSvc.h" // Event Model related classes -#include "eicd/CalorimeterHitCollection.h" -#include "eicd/RawCalorimeterHitCollection.h" +#include "edm4eic/CalorimeterHitCollection.h" +#include "edm4eic/RawCalorimeterHitCollection.h" using namespace Gaudi::Units; @@ -60,9 +60,9 @@ class ImagingPixelReco : public GaudiAlgorithm { double dyRangeADC{0}; // hits containers - DataHandle m_inputHitCollection{"inputHitCollection", Gaudi::DataHandle::Reader, + DataHandle m_inputHitCollection{"inputHitCollection", Gaudi::DataHandle::Reader, this}; - DataHandle m_outputHitCollection{"outputHitCollection", Gaudi::DataHandle::Writer, + DataHandle m_outputHitCollection{"outputHitCollection", Gaudi::DataHandle::Writer, this}; // Pointer to the geometry service @@ -145,14 +145,14 @@ class ImagingPixelReco : public GaudiAlgorithm { // create const vectors for passing to hit initializer list - const decltype(eicd::CalorimeterHitData::position) position( + const decltype(edm4eic::CalorimeterHitData::position) position( gpos.x() / m_lUnit, gpos.y() / m_lUnit, gpos.z() / m_lUnit ); - const decltype(eicd::CalorimeterHitData::local) local( + const decltype(edm4eic::CalorimeterHitData::local) local( pos.x() / m_lUnit, pos.y() / m_lUnit, pos.z() / m_lUnit ); - hits.push_back(eicd::CalorimeterHit{id, // cellID + hits.push_back(edm4eic::CalorimeterHit{id, // cellID static_cast(energy), // energy 0, // energyError static_cast(time), // time diff --git a/JugReco/src/components/ImagingTopoCluster.cpp b/JugReco/src/components/ImagingTopoCluster.cpp index 39cd0f7..e11781a 100644 --- a/JugReco/src/components/ImagingTopoCluster.cpp +++ b/JugReco/src/components/ImagingTopoCluster.cpp @@ -30,9 +30,9 @@ #include "JugReco/ClusterTypes.h" // Event Model related classes -#include "eicd/CalorimeterHitCollection.h" -#include "eicd/ProtoClusterCollection.h" -#include "eicd/vector_utils.h" +#include "edm4eic/CalorimeterHitCollection.h" +#include "edm4eic/ProtoClusterCollection.h" +#include "edm4eic/vector_utils.h" using namespace Gaudi::Units; @@ -68,10 +68,10 @@ class ImagingTopoCluster : public GaudiAlgorithm { // minimum number of hits (to save this cluster) Gaudi::Property m_minClusterNhits{this, "minClusterNhits", 10}; // input hits collection - DataHandle m_inputHitCollection{"inputHitCollection", Gaudi::DataHandle::Reader, + DataHandle m_inputHitCollection{"inputHitCollection", Gaudi::DataHandle::Reader, this}; // output clustered hits - DataHandle m_outputProtoClusterCollection{"outputProtoClusterCollection", + DataHandle m_outputProtoClusterCollection{"outputProtoClusterCollection", Gaudi::DataHandle::Writer, this}; // unitless counterparts of the input parameters @@ -135,7 +135,7 @@ class ImagingTopoCluster : public GaudiAlgorithm { // group neighboring hits std::vector visits(hits.size(), false); - std::vector>> groups; + std::vector>> groups; for (size_t i = 0; i < hits.size(); ++i) { if (msgLevel(MSG::DEBUG)) { debug() << fmt::format("hit {:d}: local position = ({}, {}, {}), global position = ({}, {}, {})", i + 1, @@ -184,7 +184,7 @@ class ImagingTopoCluster : public GaudiAlgorithm { template static inline T pow2(const T& x) { return x * x; } // helper function to group hits - bool is_neighbor(const eicd::CalorimeterHit& h1, const eicd::CalorimeterHit& h2) const { + bool is_neighbor(const edm4eic::CalorimeterHit& h1, const edm4eic::CalorimeterHit& h2) const { // different sectors, simple distance check if (h1.getSector() != h2.getSector()) { return std::sqrt(pow2(h1.getPosition().x - h2.getPosition().x) + pow2(h1.getPosition().y - h2.getPosition().y) + @@ -198,8 +198,8 @@ class ImagingTopoCluster : public GaudiAlgorithm { return (std::abs(h1.getLocal().x - h2.getLocal().x) <= localDistXY[0]) && (std::abs(h1.getLocal().y - h2.getLocal().y) <= localDistXY[1]); } else if (ldiff <= m_neighbourLayersRange) { - return (std::abs(eicd::eta(h1.getPosition()) - eicd::eta(h2.getPosition())) <= layerDistEtaPhi[0]) && - (std::abs(eicd::angleAzimuthal(h1.getPosition()) - eicd::angleAzimuthal(h2.getPosition())) <= + return (std::abs(edm4eic::eta(h1.getPosition()) - edm4eic::eta(h2.getPosition())) <= layerDistEtaPhi[0]) && + (std::abs(edm4eic::angleAzimuthal(h1.getPosition()) - edm4eic::angleAzimuthal(h2.getPosition())) <= layerDistEtaPhi[1]); } @@ -208,8 +208,8 @@ class ImagingTopoCluster : public GaudiAlgorithm { } // grouping function with Depth-First Search - void dfs_group(std::vector>& group, int idx, - const eicd::CalorimeterHitCollection& hits, std::vector& visits) const { + void dfs_group(std::vector>& group, int idx, + const edm4eic::CalorimeterHitCollection& hits, std::vector& visits) const { // not a qualified hit to participate in clustering, stop here if (hits[idx].getEnergy() < minClusterHitEdep) { visits[idx] = true; diff --git a/JugReco/src/components/InclusiveKinematicsDA.cpp b/JugReco/src/components/InclusiveKinematicsDA.cpp index fb165d6..0128656 100644 --- a/JugReco/src/components/InclusiveKinematicsDA.cpp +++ b/JugReco/src/components/InclusiveKinematicsDA.cpp @@ -21,9 +21,9 @@ using ROOT::Math::PxPyPzEVector; // Event Model related classes #include "edm4hep/MCParticleCollection.h" -#include "eicd/MCRecoParticleAssociationCollection.h" -#include "eicd/ReconstructedParticleCollection.h" -#include "eicd/InclusiveKinematicsCollection.h" +#include "edm4eic/MCRecoParticleAssociationCollection.h" +#include "edm4eic/ReconstructedParticleCollection.h" +#include "edm4eic/InclusiveKinematicsCollection.h" namespace Jug::Reco { @@ -33,15 +33,15 @@ class InclusiveKinematicsDA : public GaudiAlgorithm { "inputMCParticles", Gaudi::DataHandle::Reader, this}; - DataHandle m_inputParticleCollection{ + DataHandle m_inputParticleCollection{ "inputReconstructedParticles", Gaudi::DataHandle::Reader, this}; - DataHandle m_inputParticleAssociation{ + DataHandle m_inputParticleAssociation{ "inputParticleAssociations", Gaudi::DataHandle::Reader, this}; - DataHandle m_outputInclusiveKinematicsCollection{ + DataHandle m_outputInclusiveKinematicsCollection{ "outputInclusiveKinematics", Gaudi::DataHandle::Writer, this}; diff --git a/JugReco/src/components/InclusiveKinematicsElectron.cpp b/JugReco/src/components/InclusiveKinematicsElectron.cpp index e16ad98..af36f9e 100644 --- a/JugReco/src/components/InclusiveKinematicsElectron.cpp +++ b/JugReco/src/components/InclusiveKinematicsElectron.cpp @@ -21,9 +21,9 @@ using ROOT::Math::PxPyPzEVector; // Event Model related classes #include "edm4hep/MCParticleCollection.h" -#include "eicd/MCRecoParticleAssociationCollection.h" -#include "eicd/ReconstructedParticleCollection.h" -#include "eicd/InclusiveKinematicsCollection.h" +#include "edm4eic/MCRecoParticleAssociationCollection.h" +#include "edm4eic/ReconstructedParticleCollection.h" +#include "edm4eic/InclusiveKinematicsCollection.h" namespace Jug::Reco { @@ -33,15 +33,15 @@ class InclusiveKinematicsElectron : public GaudiAlgorithm { "inputMCParticles", Gaudi::DataHandle::Reader, this}; - DataHandle m_inputParticleCollection{ + DataHandle m_inputParticleCollection{ "inputReconstructedParticles", Gaudi::DataHandle::Reader, this}; - DataHandle m_inputParticleAssociation{ + DataHandle m_inputParticleAssociation{ "inputParticleAssociations", Gaudi::DataHandle::Reader, this}; - DataHandle m_outputInclusiveKinematicsCollection{ + DataHandle m_outputInclusiveKinematicsCollection{ "outputInclusiveKinematics", Gaudi::DataHandle::Writer, this}; diff --git a/JugReco/src/components/InclusiveKinematicsJB.cpp b/JugReco/src/components/InclusiveKinematicsJB.cpp index 3eb9af3..d4dd4cf 100644 --- a/JugReco/src/components/InclusiveKinematicsJB.cpp +++ b/JugReco/src/components/InclusiveKinematicsJB.cpp @@ -21,9 +21,9 @@ using ROOT::Math::PxPyPzEVector; // Event Model related classes #include "edm4hep/MCParticleCollection.h" -#include "eicd/MCRecoParticleAssociationCollection.h" -#include "eicd/ReconstructedParticleCollection.h" -#include "eicd/InclusiveKinematicsCollection.h" +#include "edm4eic/MCRecoParticleAssociationCollection.h" +#include "edm4eic/ReconstructedParticleCollection.h" +#include "edm4eic/InclusiveKinematicsCollection.h" namespace Jug::Reco { @@ -33,15 +33,15 @@ class InclusiveKinematicsJB : public GaudiAlgorithm { "inputMCParticles", Gaudi::DataHandle::Reader, this}; - DataHandle m_inputParticleCollection{ + DataHandle m_inputParticleCollection{ "inputReconstructedParticles", Gaudi::DataHandle::Reader, this}; - DataHandle m_inputParticleAssociation{ + DataHandle m_inputParticleAssociation{ "inputParticleAssociations", Gaudi::DataHandle::Reader, this}; - DataHandle m_outputInclusiveKinematicsCollection{ + DataHandle m_outputInclusiveKinematicsCollection{ "outputInclusiveKinematics", Gaudi::DataHandle::Writer, this}; diff --git a/JugReco/src/components/InclusiveKinematicsSigma.cpp b/JugReco/src/components/InclusiveKinematicsSigma.cpp index 8070171..03c2156 100644 --- a/JugReco/src/components/InclusiveKinematicsSigma.cpp +++ b/JugReco/src/components/InclusiveKinematicsSigma.cpp @@ -21,9 +21,9 @@ using ROOT::Math::PxPyPzEVector; // Event Model related classes #include "edm4hep/MCParticleCollection.h" -#include "eicd/MCRecoParticleAssociationCollection.h" -#include "eicd/ReconstructedParticleCollection.h" -#include "eicd/InclusiveKinematicsCollection.h" +#include "edm4eic/MCRecoParticleAssociationCollection.h" +#include "edm4eic/ReconstructedParticleCollection.h" +#include "edm4eic/InclusiveKinematicsCollection.h" namespace Jug::Reco { @@ -33,15 +33,15 @@ class InclusiveKinematicsSigma : public GaudiAlgorithm { "inputMCParticles", Gaudi::DataHandle::Reader, this}; - DataHandle m_inputParticleCollection{ + DataHandle m_inputParticleCollection{ "inputReconstructedParticles", Gaudi::DataHandle::Reader, this}; - DataHandle m_inputParticleAssociation{ + DataHandle m_inputParticleAssociation{ "inputParticleAssociations", Gaudi::DataHandle::Reader, this}; - DataHandle m_outputInclusiveKinematicsCollection{ + DataHandle m_outputInclusiveKinematicsCollection{ "outputInclusiveKinematics", Gaudi::DataHandle::Writer, this}; diff --git a/JugReco/src/components/InclusiveKinematicseSigma.cpp b/JugReco/src/components/InclusiveKinematicseSigma.cpp index 1fca536..f9c5285 100644 --- a/JugReco/src/components/InclusiveKinematicseSigma.cpp +++ b/JugReco/src/components/InclusiveKinematicseSigma.cpp @@ -21,9 +21,9 @@ using ROOT::Math::PxPyPzEVector; // Event Model related classes #include "edm4hep/MCParticleCollection.h" -#include "eicd/MCRecoParticleAssociationCollection.h" -#include "eicd/ReconstructedParticleCollection.h" -#include "eicd/InclusiveKinematicsCollection.h" +#include "edm4eic/MCRecoParticleAssociationCollection.h" +#include "edm4eic/ReconstructedParticleCollection.h" +#include "edm4eic/InclusiveKinematicsCollection.h" namespace Jug::Reco { @@ -33,15 +33,15 @@ class InclusiveKinematicseSigma : public GaudiAlgorithm { "inputMCParticles", Gaudi::DataHandle::Reader, this}; - DataHandle m_inputParticleCollection{ + DataHandle m_inputParticleCollection{ "inputReconstructedParticles", Gaudi::DataHandle::Reader, this}; - DataHandle m_inputParticleAssociation{ + DataHandle m_inputParticleAssociation{ "inputParticleAssociations", Gaudi::DataHandle::Reader, this}; - DataHandle m_outputInclusiveKinematicsCollection{ + DataHandle m_outputInclusiveKinematicsCollection{ "outputInclusiveKinematics", Gaudi::DataHandle::Writer, this}; diff --git a/JugReco/src/components/ParticleCollector.cpp b/JugReco/src/components/ParticleCollector.cpp index 24e28b5..59f4306 100644 --- a/JugReco/src/components/ParticleCollector.cpp +++ b/JugReco/src/components/ParticleCollector.cpp @@ -10,7 +10,7 @@ #include "JugBase/DataHandle.h" // Event Model related classes -#include "eicd/ReconstructedParticleCollection.h" +#include "edm4eic/ReconstructedParticleCollection.h" namespace Jug::Reco { @@ -24,10 +24,10 @@ namespace Jug::Reco { class ParticleCollector : public GaudiAlgorithm { private: Gaudi::Property> m_inputParticles{this, "inputParticles", {}, "Particles to be aggregated"}; - DataHandle m_outputParticles{"outputParticles", Gaudi::DataHandle::Writer, + DataHandle m_outputParticles{"outputParticles", Gaudi::DataHandle::Writer, this}; - std::vector*> m_particleCollections; + std::vector*> m_particleCollections; public: ParticleCollector(const std::string& name, ISvcLocator* svcLoc) : GaudiAlgorithm(name, svcLoc) { @@ -46,7 +46,7 @@ class ParticleCollector : public GaudiAlgorithm { for (auto colname : m_inputParticles) { debug() << "initializing collection: " << colname << endmsg; m_particleCollections.push_back( - new DataHandle{colname, Gaudi::DataHandle::Reader, this}); + new DataHandle{colname, Gaudi::DataHandle::Reader, this}); } return StatusCode::SUCCESS; } diff --git a/JugReco/src/components/PhotoMultiplierReco.cpp b/JugReco/src/components/PhotoMultiplierReco.cpp index 141446b..652a044 100644 --- a/JugReco/src/components/PhotoMultiplierReco.cpp +++ b/JugReco/src/components/PhotoMultiplierReco.cpp @@ -28,8 +28,8 @@ #include "JugBase/IGeoSvc.h" // Event Model related classes -#include "eicd/PMTHitCollection.h" -#include "eicd/RawPMTHitCollection.h" +#include "edm4eic/PMTHitCollection.h" +#include "edm4eic/RawPMTHitCollection.h" using namespace Gaudi::Units; @@ -44,8 +44,8 @@ namespace Jug::Reco { */ class PhotoMultiplierReco : public GaudiAlgorithm { private: - DataHandle m_inputHitCollection{"inputHitCollection", Gaudi::DataHandle::Reader, this}; - DataHandle m_outputHitCollection{"outputHitCollection", Gaudi::DataHandle::Writer, this}; + DataHandle m_inputHitCollection{"inputHitCollection", Gaudi::DataHandle::Reader, this}; + DataHandle m_outputHitCollection{"outputHitCollection", Gaudi::DataHandle::Writer, this}; Gaudi::Property m_timeStep{this, "timeStep", 0.0625 * ns}; Gaudi::Property m_minNpe{this, "minNpe", 0.0}; Gaudi::Property m_speMean{this, "speMean", 80.0}; @@ -91,7 +91,7 @@ class PhotoMultiplierReco : public GaudiAlgorithm { auto pos = m_geoSvc->cellIDPositionConverter()->findContext(id)->volumePlacement().position(); // cell dimension auto dim = m_geoSvc->cellIDPositionConverter()->cellDimensions(id); - hits.push_back(eicd::PMTHit{ + hits.push_back(edm4eic::PMTHit{ rh.getCellID(), npe, time, diff --git a/JugReco/src/components/SimpleClustering.cpp b/JugReco/src/components/SimpleClustering.cpp index 446a736..a48b424 100644 --- a/JugReco/src/components/SimpleClustering.cpp +++ b/JugReco/src/components/SimpleClustering.cpp @@ -20,11 +20,11 @@ // Event Model related classes #include "edm4hep/SimCalorimeterHitCollection.h" -#include "eicd/CalorimeterHitCollection.h" -#include "eicd/ClusterCollection.h" -#include "eicd/ProtoClusterCollection.h" -#include "eicd/RawCalorimeterHitCollection.h" -#include "eicd/vector_utils.h" +#include "edm4eic/CalorimeterHitCollection.h" +#include "edm4eic/ClusterCollection.h" +#include "edm4eic/ProtoClusterCollection.h" +#include "edm4eic/RawCalorimeterHitCollection.h" +#include "edm4eic/vector_utils.h" using namespace Gaudi::Units; @@ -36,9 +36,9 @@ namespace Jug::Reco { */ class SimpleClustering : public GaudiAlgorithm { private: - using RecHits = eicd::CalorimeterHitCollection; - using ProtoClusters = eicd::ProtoClusterCollection; - using Clusters = eicd::ClusterCollection; + using RecHits = edm4eic::CalorimeterHitCollection; + using ProtoClusters = edm4eic::ProtoClusterCollection; + using Clusters = edm4eic::ClusterCollection; DataHandle m_inputHitCollection{"inputHitCollection", Gaudi::DataHandle::Reader, this}; DataHandle m_outputProtoClusters{"outputProtoCluster", Gaudi::DataHandle::Writer, this}; @@ -97,13 +97,13 @@ namespace Jug::Reco { // mcHits = m_inputMC->get(); //} - std::vector> the_hits; - std::vector> remaining_hits; + std::vector> the_hits; + std::vector> remaining_hits; double max_dist = m_maxDistance.value() / mm; double min_energy = m_minModuleEdep.value() / GeV; - eicd::CalorimeterHit ref_hit; + edm4eic::CalorimeterHit ref_hit; bool have_ref = false; // Collect all our hits, and get the highest energy hit { @@ -124,10 +124,10 @@ namespace Jug::Reco { while (have_ref && ref_hit.getEnergy() > min_energy) { - std::vector> cluster_hits; + std::vector> cluster_hits; for (const auto& [idx, h] : the_hits) { - if (eicd::magnitude(h.getPosition() - ref_hit.getPosition()) < max_dist) { + if (edm4eic::magnitude(h.getPosition() - ref_hit.getPosition()) < max_dist) { cluster_hits.emplace_back(idx, h); } else { remaining_hits.emplace_back(idx, h); @@ -136,7 +136,7 @@ namespace Jug::Reco { double total_energy = std::accumulate( std::begin(cluster_hits), std::end(cluster_hits), 0.0, - [](double t, const std::pair& h1) { return (t + h1.second.getEnergy()); }); + [](double t, const std::pair& h1) { return (t + h1.second.getEnergy()); }); if (msgLevel(MSG::DEBUG)) { debug() << " total_energy = " << total_energy << endmsg; diff --git a/JugReco/src/components/TrackerHitReconstruction.cpp b/JugReco/src/components/TrackerHitReconstruction.cpp index b669208..d054f7b 100644 --- a/JugReco/src/components/TrackerHitReconstruction.cpp +++ b/JugReco/src/components/TrackerHitReconstruction.cpp @@ -21,8 +21,8 @@ // Event Model related classes //#include "GaudiExamples/MyTrack.h" -#include "eicd/RawTrackerHitCollection.h" -#include "eicd/TrackerHitCollection.h" +#include "edm4eic/RawTrackerHitCollection.h" +#include "edm4eic/TrackerHitCollection.h" #include "DD4hep/DD4hepUnits.h" @@ -46,9 +46,9 @@ namespace Jug::Reco { class TrackerHitReconstruction : public GaudiAlgorithm { private: Gaudi::Property m_timeResolution{this, "timeResolution", 10}; // in ns - DataHandle m_inputHitCollection{"inputHitCollection", Gaudi::DataHandle::Reader, + DataHandle m_inputHitCollection{"inputHitCollection", Gaudi::DataHandle::Reader, this}; - DataHandle m_outputHitCollection{"outputHitCollection", Gaudi::DataHandle::Writer, + DataHandle m_outputHitCollection{"outputHitCollection", Gaudi::DataHandle::Writer, this}; /// Pointer to the geometry service @@ -106,7 +106,7 @@ namespace Jug::Reco { // - XYZ segmentation: xx -> sigma_x, yy-> sigma_y, zz -> sigma_z, tt -> 0 // This is properly in line with how we get the local coordinates for the hit // in the TrackerSourceLinker. - eicd::TrackerHit hit{ahit.getCellID(), // Raw DD4hep cell ID + edm4eic::TrackerHit hit{ahit.getCellID(), // Raw DD4hep cell ID {static_cast(pos.x() / mm), static_cast(pos.y() / mm), static_cast(pos.z() / mm)}, // mm {get_variance(dim[0] / mm), get_variance(dim[1] / mm), // variance (see note above) diff --git a/JugReco/src/components/TrackingHitsCollector.cpp b/JugReco/src/components/TrackingHitsCollector.cpp index 2172314..0e4d1ea 100644 --- a/JugReco/src/components/TrackingHitsCollector.cpp +++ b/JugReco/src/components/TrackingHitsCollector.cpp @@ -10,7 +10,7 @@ #include "JugBase/DataHandle.h" // Event Model related classes -#include "eicd/TrackerHitCollection.h" +#include "edm4eic/TrackerHitCollection.h" namespace Jug::Reco { @@ -20,12 +20,12 @@ namespace Jug::Reco { */ class TrackingHitsCollector : public GaudiAlgorithm { private: - DataHandle m_trackerBarrelHits{"trackerBarrelHits", Gaudi::DataHandle::Reader, this}; - DataHandle m_trackerEndcapHits{"trackerEndcapHits", Gaudi::DataHandle::Reader, this}; - DataHandle m_vertexBarrelHits {"vertexBarrelHits" , Gaudi::DataHandle::Reader, this}; - DataHandle m_vertexEndcapHits {"vertexEndcapHits" , Gaudi::DataHandle::Reader, this}; - DataHandle m_gemEndcapHits {"gemEndcapHits" , Gaudi::DataHandle::Reader, this}; - DataHandle m_outputHitCollection{"outputHitCollection", Gaudi::DataHandle::Writer, this}; + DataHandle m_trackerBarrelHits{"trackerBarrelHits", Gaudi::DataHandle::Reader, this}; + DataHandle m_trackerEndcapHits{"trackerEndcapHits", Gaudi::DataHandle::Reader, this}; + DataHandle m_vertexBarrelHits {"vertexBarrelHits" , Gaudi::DataHandle::Reader, this}; + DataHandle m_vertexEndcapHits {"vertexEndcapHits" , Gaudi::DataHandle::Reader, this}; + DataHandle m_gemEndcapHits {"gemEndcapHits" , Gaudi::DataHandle::Reader, this}; + DataHandle m_outputHitCollection{"outputHitCollection", Gaudi::DataHandle::Writer, this}; public: TrackingHitsCollector(const std::string& name, ISvcLocator* svcLoc) @@ -47,11 +47,11 @@ namespace Jug::Reco { StatusCode execute() override { - const eicd::TrackerHitCollection* trkBarrelHits = m_trackerBarrelHits.get(); - const eicd::TrackerHitCollection* trkEndcapHits = m_trackerEndcapHits.get(); - const eicd::TrackerHitCollection* vtxBarrelHits = m_vertexBarrelHits .get(); - const eicd::TrackerHitCollection* vtxEndcapHits = m_vertexEndcapHits .get(); - const eicd::TrackerHitCollection* gemEndcapHits = m_gemEndcapHits .get(); + const edm4eic::TrackerHitCollection* trkBarrelHits = m_trackerBarrelHits.get(); + const edm4eic::TrackerHitCollection* trkEndcapHits = m_trackerEndcapHits.get(); + const edm4eic::TrackerHitCollection* vtxBarrelHits = m_vertexBarrelHits .get(); + const edm4eic::TrackerHitCollection* vtxEndcapHits = m_vertexEndcapHits .get(); + const edm4eic::TrackerHitCollection* gemEndcapHits = m_gemEndcapHits .get(); auto* outputHits = m_outputHitCollection.createAndPut(); for (const auto* hits : {trkBarrelHits, trkEndcapHits, vtxBarrelHits, vtxEndcapHits, gemEndcapHits}) { diff --git a/JugReco/src/components/TrackingHitsCollector2.cpp b/JugReco/src/components/TrackingHitsCollector2.cpp index da51534..17afb4f 100644 --- a/JugReco/src/components/TrackingHitsCollector2.cpp +++ b/JugReco/src/components/TrackingHitsCollector2.cpp @@ -10,7 +10,7 @@ #include "JugBase/DataHandle.h" // Event Model related classes -#include "eicd/TrackerHitCollection.h" +#include "edm4eic/TrackerHitCollection.h" namespace Jug::Reco { @@ -24,9 +24,9 @@ namespace Jug::Reco { class TrackingHitsCollector2 : public GaudiAlgorithm { private: Gaudi::Property> m_inputTrackingHits{this, "inputTrackingHits", {},"Tracker hits to be aggregated"}; - DataHandle m_trackingHits{"trackingHits", Gaudi::DataHandle::Writer, this}; + DataHandle m_trackingHits{"trackingHits", Gaudi::DataHandle::Writer, this}; - std::vector*> m_hitCollections; + std::vector*> m_hitCollections; public: TrackingHitsCollector2(const std::string& name, ISvcLocator* svcLoc) @@ -46,7 +46,7 @@ namespace Jug::Reco { } for (auto colname : m_inputTrackingHits) { debug() << "initializing collection: " << colname << endmsg; - m_hitCollections.push_back(new DataHandle{colname, Gaudi::DataHandle::Reader, this}); + m_hitCollections.push_back(new DataHandle{colname, Gaudi::DataHandle::Reader, this}); } return StatusCode::SUCCESS; } @@ -58,7 +58,7 @@ namespace Jug::Reco { debug() << "execute collector" << endmsg; } for(const auto& hits: m_hitCollections) { - const eicd::TrackerHitCollection* hitCol = hits->get(); + const edm4eic::TrackerHitCollection* hitCol = hits->get(); if (msgLevel(MSG::DEBUG)) { debug() << "col n hits: " << hitCol->size() << endmsg; } diff --git a/JugTrack/CMakeLists.txt b/JugTrack/CMakeLists.txt index fc6bf33..47a3beb 100644 --- a/JugTrack/CMakeLists.txt +++ b/JugTrack/CMakeLists.txt @@ -14,17 +14,15 @@ gaudi_add_module(JugTrackPlugins JugBase ROOT::Core ROOT::RIO ROOT::Tree EDM4HEP::edm4hep - EICD::eicd + EDM4EIC::edm4eic DD4hep::DDRec ActsCore - ${genfit2} ) target_include_directories(JugTrackPlugins PUBLIC $ $ $ - ${genfit2_INCLUDE_DIR} ) target_compile_options(JugTrackPlugins PRIVATE -Wno-suggest-override) diff --git a/JugTrack/JugTrack/SourceLinks.h b/JugTrack/JugTrack/SourceLinks.h index 898ace3..0598eda 100644 --- a/JugTrack/JugTrack/SourceLinks.h +++ b/JugTrack/JugTrack/SourceLinks.h @@ -10,7 +10,7 @@ #include #include -#include "eicd/TrackerHitCollection.h" +#include "edm4eic/TrackerHitCollection.h" namespace Jug { @@ -34,7 +34,7 @@ class SourceLink { // need to store pointers to make the object copyable const Acts::Surface* m_surface; //const ActsFatras::Hit* m_truthHit; - const eicd::TrackerHit* m_Hit ; + const edm4eic::TrackerHit* m_Hit ; public: SourceLink(const Acts::Surface& surface, //const ActsFatras::Hit& truthHit, diff --git a/JugTrack/src/components/CKFTracking.cpp b/JugTrack/src/components/CKFTracking.cpp index adabc3d..2875edc 100644 --- a/JugTrack/src/components/CKFTracking.cpp +++ b/JugTrack/src/components/CKFTracking.cpp @@ -40,7 +40,7 @@ #include "JugTrack/Track.hpp" -#include "eicd/TrackerHitCollection.h" +#include "edm4eic/TrackerHitCollection.h" #include #include diff --git a/JugTrack/src/components/CKFTracking.h b/JugTrack/src/components/CKFTracking.h index ead37d2..aff0380 100644 --- a/JugTrack/src/components/CKFTracking.h +++ b/JugTrack/src/components/CKFTracking.h @@ -23,7 +23,7 @@ #include "JugTrack/Track.hpp" #include "JugTrack/Trajectories.hpp" -#include "eicd/TrackerHitCollection.h" +#include "edm4eic/TrackerHitCollection.h" #include "Acts/Definitions/Common.hpp" #include "Acts/Geometry/TrackingGeometry.hpp" diff --git a/JugTrack/src/components/ConformalXYPeakProtoTracks.cpp b/JugTrack/src/components/ConformalXYPeakProtoTracks.cpp index 5ee7719..94a9da2 100644 --- a/JugTrack/src/components/ConformalXYPeakProtoTracks.cpp +++ b/JugTrack/src/components/ConformalXYPeakProtoTracks.cpp @@ -18,7 +18,7 @@ #include "Math/Vector3D.h" #include "Math/Vector2D.h" -#include "eicd/TrackerHitCollection.h" +#include "edm4eic/TrackerHitCollection.h" namespace Jug::Reco { @@ -30,7 +30,7 @@ namespace Jug::Reco { */ class ConformalXYPeakProtoTracks : public GaudiAlgorithm { private: - DataHandle m_inputTrackerHits{"inputTrackerHits", Gaudi::DataHandle::Reader, this}; + DataHandle m_inputTrackerHits{"inputTrackerHits", Gaudi::DataHandle::Reader, this}; DataHandle m_outputProtoTracks{"outputProtoTracks", Gaudi::DataHandle::Writer, this}; DataHandle m_nProtoTracks{"nProtoTracks", Gaudi::DataHandle::Writer, this}; @@ -54,7 +54,7 @@ class ConformalXYPeakProtoTracks : public GaudiAlgorithm { StatusCode execute() override { // input collection - const eicd::TrackerHitCollection* hits = m_inputTrackerHits.get(); + const edm4eic::TrackerHitCollection* hits = m_inputTrackerHits.get(); // Create output collections auto* proto_tracks = m_outputProtoTracks.createAndPut(); int n_proto_tracks = 0; diff --git a/JugTrack/src/components/FinderAlgoTemplate.cpp b/JugTrack/src/components/FinderAlgoTemplate.cpp index 0425abd..c2238d4 100644 --- a/JugTrack/src/components/FinderAlgoTemplate.cpp +++ b/JugTrack/src/components/FinderAlgoTemplate.cpp @@ -16,7 +16,7 @@ #include "Math/Vector3D.h" -#include "eicd/TrackerHitCollection.h" +#include "edm4eic/TrackerHitCollection.h" namespace Jug::Reco { @@ -26,7 +26,7 @@ namespace Jug::Reco { */ class FinderAlgoTemplate : public GaudiAlgorithm { private: - DataHandle m_inputTrackerHits{"inputTrackerHits", Gaudi::DataHandle::Reader, this}; + DataHandle m_inputTrackerHits{"inputTrackerHits", Gaudi::DataHandle::Reader, this}; DataHandle m_outputProtoTracks{"outputProtoTracks", Gaudi::DataHandle::Writer, this}; public: @@ -44,7 +44,7 @@ class FinderAlgoTemplate : public GaudiAlgorithm { StatusCode execute() override { // input collection - //const eicd::TrackerHitCollection* hits = m_inputTrackerHits.get(); + //const edm4eic::TrackerHitCollection* hits = m_inputTrackerHits.get(); // Create output collections /*auto proto_tracks = */m_outputProtoTracks.createAndPut(); diff --git a/JugTrack/src/components/GenFitTrackFitter.cpp b/JugTrack/src/components/GenFitTrackFitter.cpp deleted file mode 100644 index f8c52aa..0000000 --- a/JugTrack/src/components/GenFitTrackFitter.cpp +++ /dev/null @@ -1,308 +0,0 @@ -// SPDX-License-Identifier: LGPL-3.0-or-later -// Copyright (C) 2022 Whitney Armstrong, Sylvester Joosten - -#include "GenFitTrackFitter.h" -// Gaudi -#include "Gaudi/Property.h" -#include "GaudiAlg/GaudiAlgorithm.h" -#include "GaudiAlg/GaudiTool.h" -#include "GaudiAlg/Transformer.h" -#include "GaudiKernel/RndmGenerators.h" -#include "GaudiKernel/ToolHandle.h" - -#include "DDRec/CellIDPositionConverter.h" -#include "DDRec/Surface.h" -#include "DDRec/SurfaceManager.h" - -#include "JugBase/BField/DD4hepBField.h" -#include "JugBase/DataHandle.h" -#include "JugBase/IGeoSvc.h" -#include "JugTrack/GeometryContainers.hpp" -#include "JugTrack/IndexSourceLink.hpp" -#include "JugTrack/Measurement.hpp" -#include "JugTrack/Track.hpp" - -#include "eicd/TrackerHitCollection.h" -#include "eicd/vector_utils.h" - -#include -#include -#include -#include - -//# genfit -#include "ConstField.h" -#include "DAF.h" -#include "Exception.h" -#include "FieldManager.h" -#include "KalmanFitterRefTrack.h" -#include "MaterialEffects.h" -#include "RKTrackRep.h" -#include "StateOnPlane.h" -#include "TGeoMaterialInterface.h" -#include "Track.h" -#include "TrackPoint.h" -//#include -#include "HelixTrackModel.h" -#include "PlanarMeasurement.h" -//#include "MeasurementCreator.h" -#include "WireMeasurement.h" - -#include "TDatabasePDG.h" -#include "TEveManager.h" -#include "TGeoManager.h" -#include "TRandom.h" -#include "TVector3.h" -#include - -namespace Jug::Reco { - -GenFitTrackFitter::GenFitTrackFitter(const std::string& name, ISvcLocator* svcLoc) : GaudiAlgorithm(name, svcLoc) { - declareProperty("inputHitCollection", m_inputHitCollection, ""); - declareProperty("initialTrackParameters", m_initialTrackParameters, ""); - declareProperty("inputProtoTracks", m_inputProtoTracks, ""); - declareProperty("trackParameters", m_foundTracks, ""); - declareProperty("outputTrajectories", m_outputTrajectories, ""); -} - -StatusCode GenFitTrackFitter::initialize() { - if (GaudiAlgorithm::initialize().isFailure()) { - return StatusCode::FAILURE; - } - 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; - } - - genfit::FieldManager::getInstance()->init(new FieldImp(m_geoSvc->detector())); - // 0., 0., m_geoSvc->centralMagneticField() * 10.0)); // gentfit uses kilo-Gauss - genfit::MaterialEffects::getInstance()->init(new genfit::TGeoMaterialInterface()); - - // copy the whole map to get around genfit's interface - // this should be returning a const& - m_detPlaneMap = m_geoSvc->getDetPlaneMap(); - m_surfaceMap = m_geoSvc->getDD4hepSurfaceMap(); - - return StatusCode::SUCCESS; -} - -StatusCode GenFitTrackFitter::execute() { - // Read input data - const eicd::TrackerHitCollection* hits = m_inputHitCollection.get(); - const TrackParametersContainer* initialParameters = m_initialTrackParameters.get(); - const ProtoTrackContainer* protoTracks = m_inputProtoTracks.get(); - - // TrajectoryContainer trajectories; - // commented out unused variables - /*auto trajectories =*/m_outputTrajectories.createAndPut(); - /*auto trackParameters =*/m_foundTracks.createAndPut(); - - int n_tracks = initialParameters->size(); - int n_proto_tracks = protoTracks->size(); - // Unused variable - // int ID = 0; - - // Assuming init track parameters have been match with proto tracks by index - if (n_proto_tracks != n_tracks) { - warning() << " Number of proto tracks does not match the initial track parameters." << endmsg; - } - - for (int itrack = 0; itrack < std::min(n_tracks, n_proto_tracks); itrack++) { - const auto& track_param = (*initialParameters)[itrack]; - const auto& proto_track = (*protoTracks)[itrack]; - - if (msgLevel(MSG::DEBUG)) { - debug() << "track mom : " << track_param.absoluteMomentum() << endmsg; - } - if (hits->size() < 2) { - return StatusCode::SUCCESS; - } - // init fitter - // genfit::KalmanFitterRefTrack fitter; - // fitter.setDebugLvl(1); - genfit::DAF fitter; - // genfit::AbsKalmanFitter* fitter = new genfit::KalmanFitterRefTrack(); - - ROOT::Math::XYZVector tp(track_param.momentum()[0], track_param.momentum()[1], track_param.momentum()[2]); - auto first_hit = (*hits)[proto_track[0]]; - auto first_hit_phi = eicd::angleAzimuthal(first_hit.getPosition()); - auto track_param_phi = tp.phi(); - if (msgLevel(MSG::DEBUG)) { - debug() << " first hit phi: " << first_hit_phi << endmsg; - debug() << "init track phi: " << track_param_phi << endmsg; - } - if (std::fabs(first_hit_phi - track_param_phi) > 0.15) { - warning() << "Seed directions does not match first hit phi. " << endmsg; - continue; - } - - // start values for the fit, e.g. from pattern recognition - TVector3 pos(0, 0, 0); - TVector3 mom(track_param.momentum()[0], track_param.momentum()[1], track_param.momentum()[2]); - TMatrixDSym covM(6); - covM(0, 0) = 0.001; - covM(1, 1) = 0.001; - covM(2, 2) = 1.0; - covM(3, 3) = 0.05 * track_param.momentum()[0] * 0.05 * track_param.momentum()[0]; - covM(4, 4) = 0.05 * track_param.momentum()[1] * 0.05 * track_param.momentum()[1]; - covM(5, 5) = 0.05 * track_param.momentum()[2] * 0.05 * track_param.momentum()[2]; - - if (msgLevel(MSG::DEBUG)) { - debug() << "covM = " << covM(0, 0) << "," << covM(1, 1) << "," << covM(2, 2) << "," << covM(3, 3) << "," - << covM(4, 4) << "," << covM(5, 5) << " " << endmsg; - } - - // trackrep - // @FIXME: raw new should be avoided, either place on the stack or use - // std::unique_ptr<> - genfit::AbsTrackRep* electron_rep = new genfit::RKTrackRep(11); - // unusud - // genfit::AbsTrackRep* positron_rep = new genfit::RKTrackRep(-11); - // genfit::AbsTrackRep* piplus_rep = new genfit::RKTrackRep(211); - // genfit::AbsTrackRep* piminus_rep = new genfit::RKTrackRep(-211); - - // smeared start state - genfit::MeasuredStateOnPlane stateSmeared(electron_rep); - stateSmeared.setPosMomCov(pos, mom, covM); - - // create track - TVectorD seedState(6); - TMatrixDSym seedCov(6); - stateSmeared.get6DStateCov(seedState, seedCov); - // genfit::Track fitTrack(rep, seedState, seedCov); - - // create track - genfit::Track fitTrack(electron_rep, seedState, seedCov); - // genfit::Track fitTrack(electron_rep, pos, mom); - // fitTrack.addTrackRep(positron_rep); - // fitTrack.addTrackRep(piplus_rep ); - // fitTrack.addTrackRep(piminus_rep); - - if (msgLevel(MSG::DEBUG)) { - debug() << (*hits).size() << " hits " << endmsg; - } - - int nhit = 0; - for (int ihit : proto_track) { - const auto& ahit = (*hits)[ihit]; - - const auto* vol_ctx = m_geoSvc->cellIDPositionConverter()->findContext(ahit.getCellID()); - auto vol_id = vol_ctx->identifier; - auto volman = m_geoSvc->detector()->volumeManager(); - auto alignment = volman.lookupDetElement(vol_id).nominal(); - auto local_position = alignment.worldToLocal( - {ahit.getPosition().x / 10.0, ahit.getPosition().y / 10.0, ahit.getPosition().z / 10.0}); - auto* surf = m_surfaceMap[vol_id]; - auto local_position2 = - surf->globalToLocal({ahit.getPosition().x / 10.0, ahit.getPosition().y / 10.0, ahit.getPosition().z / 10.0}); - - TMatrixDSym hitCov(2); - hitCov.UnitMatrix(); - hitCov(0, 0) = ahit.getPositionError().xx / (100.0); // go from mm^2 to cm^2 - hitCov(1, 1) = ahit.getPositionError().yy / (100.0); // go from mm^2 to cm^2 - - if (msgLevel(MSG::DEBUG)) { - debug() << "------------------------------------ " << endmsg; - debug() << " hit position : " << ahit.getPosition().x / 10 << " " << ahit.getPosition().y / 10 << " " - << ahit.getPosition().z / 10 << endmsg; - debug() << " dd4hep loc pos : " << local_position.x() << " " << local_position.y() << " " - << local_position.z() << endmsg; - debug() << " dd4hep surf pos : " << local_position2.u() << " " << local_position2.v() << endmsg; - } - - /** \todo Add check for XZ segmentations to use the right local coordinates. - * Unlike acts, the conversion to the local system isn't going from 3D -> 2D. - * Thefore there is one coordinate that is zero. Which one depends on the the segmentation - * type XY, XZ, etc. For XY the Z-coordinate is zero. For the XZ, the Y-coordinate is zero. - */ - TVectorD hitCoords(2); - hitCoords[0] = local_position2.u(); - hitCoords[1] = local_position2.v(); - if (msgLevel(MSG::DEBUG)) { - debug() << "covariance matrix : " << hitCov(0, 0) << " " << hitCov(1, 1) << " " << endmsg; - debug() << " hit coordinates : " << hitCoords[0] << " " << hitCoords[1] << " " << endmsg; - } - auto* measurement = new genfit::PlanarMeasurement(hitCoords, hitCov, 1 /** type **/, nhit, nullptr); - - // measurement->setPlane(genfit::SharedPlanePtr(new genfit::DetPlane(point, u_dir, v_dir)), - measurement->setPlane(m_detPlaneMap[vol_id], vol_id); - fitTrack.insertPoint(new genfit::TrackPoint(measurement, &fitTrack)); - // positronFitTrack.insertPoint(new genfit::TrackPoint(measurement, &positronFitTrack)); - - nhit++; - } - - // do the fit - try { - if (msgLevel(MSG::DEBUG)) { - debug() << "Electron track: " << endmsg; - fitTrack.checkConsistency(); - } - fitter.processTrack(&fitTrack, true); - bool isConverged = fitTrack.getFitStatus()->isFitConverged(); - if (!isConverged) { - fitter.processTrack(&fitTrack, true); - } - - isConverged = fitTrack.getFitStatus()->isFitConverged(); - if (!isConverged) { - fitter.processTrack(&fitTrack); - } - - // print fit result - fitTrack.getFittedState().Print(); - // isConverged = fitTrack.getFitStatus()->isFitConverged(); - - bool isFitted = fitTrack.getFitStatus()->isFitted(); - float chi2 = fitTrack.getFitStatus()->getChi2(); - float ndf = fitTrack.getFitStatus()->getNdf(); - // unused - // float charge = fitTrack.getFitStatus()->getCharge(); - - TVector3 vertexPos; - TVector3 vertexMom; - TMatrixDSym vertexCov; - genfit::MeasuredStateOnPlane state = fitTrack.getFittedState(); // copy - TVector3 vertex(0, 0, 0); - TVector3 axis(0, 0, 1); - // state.extrapolateToPoint(vertex); - // or alternatively - state.extrapolateToLine(vertex, axis); - state.getPosMomCov(vertexPos, vertexMom, vertexCov); - - if (msgLevel(MSG::DEBUG)) { - debug() << "Electron track: " << endmsg; - fitTrack.checkConsistency(); - debug() << "vertex pos: " << vertexPos.x() << ", " << vertexPos.y() << ", " << vertexPos.z() << endmsg; - debug() << "vertex mom: " << vertexMom.x() << ", " << vertexMom.y() << ", " << vertexMom.z() << endmsg; - debug() << "track status: " << endmsg; - debug() << " fitted = " << isFitted << endmsg; - debug() << " converged =" << isConverged << endmsg; - debug() << " chi2/ndf = " << isFitted << "/" << ndf << " = " << chi2 / ndf << endmsg; - debug() << " charge =" << isConverged << endmsg; - // debug() << "Positron track: " << endmsg; - // positronFitTrack.checkConsistency(); - } - } catch (genfit::Exception& e) { - warning() << e.what() << endmsg; - warning() << "Exception, next track" << endmsg; - continue; - } - - // eicd::TrackParameters electron_track_params({ID++, algorithmID()}, {0.0,0.0},{0.0,0.0},{0.0,0.0},{0.0,0.0}, - - // TrackParameters(eicd::Index ID, eicd::FloatPair loc, eicd::FloatPair locError, eicd::Direction direction, - // eicd::Direction directionError, float qOverP, float qOverPError, float time, float timeError); - // tracks->push_back(electron_track_params); - - // delete fitter; - } - - return StatusCode::SUCCESS; -} - -// NOLINTNEXTLINE(cppcoreguidelines-avoid-non-const-global-variables) -DECLARE_COMPONENT(GenFitTrackFitter) -} // namespace Jug::Reco diff --git a/JugTrack/src/components/GenFitTrackFitter.h b/JugTrack/src/components/GenFitTrackFitter.h deleted file mode 100644 index 8f2e254..0000000 --- a/JugTrack/src/components/GenFitTrackFitter.h +++ /dev/null @@ -1,103 +0,0 @@ -// SPDX-License-Identifier: LGPL-3.0-or-later -// Copyright (C) 2022 Whitney Armstrong - -#ifndef JUGGLER_JUGRECO_GenFitTrackFitter_HH -#define JUGGLER_JUGRECO_GenFitTrackFitter_HH - -#include -#include -#include - -// Gaudi -#include "GaudiAlg/GaudiAlgorithm.h" -#include "GaudiKernel/ToolHandle.h" -#include "Gaudi/Property.h" - -#include "JugBase/DataHandle.h" -#include "JugBase/IGeoSvc.h" -#include "JugBase/BField/DD4hepBField.h" -#include "JugTrack/GeometryContainers.hpp" -#include "JugTrack/IndexSourceLink.hpp" -#include "JugTrack/Track.hpp" -#include "JugTrack/Measurement.hpp" -#include "JugTrack/Trajectories.hpp" -#include "JugTrack/ProtoTrack.hpp" - -#include "eicd/TrackerHitCollection.h" -#include "eicd/TrajectoryCollection.h" -#include "eicd/TrackParametersCollection.h" - -//genfitk -#include "FieldManager.h" - -#include -#include - -namespace Jug::Reco { - - - /** Genfit based tracking algorithm. - * - * \ingroup tracking - */ - class GenFitTrackFitter : public GaudiAlgorithm { - public: - - class FieldImp : public genfit::AbsBField { - protected: - dd4hep::Detector* m_detector; - public: - FieldImp(dd4hep::Detector* det): m_detector(det) {} - virtual ~FieldImp() {} - - /** Get the magneticField [kGauss] at position. - * - * Note that tgeo units are used. [cm] and [kGauss]. - */ - TVector3 get(const TVector3& position) const override { - double pos[3] = {position.x(), position.y(), position.z()}; - double field[3]; - this->get(pos[0], pos[1], pos[2], field[0], field[1], field[2]); - return {field[0], field[1], field[2]}; - } - - /** Get the magneticField [kGauss] at position. - * - * Note that tgeo units are used. [cm] and [kGauss]. - */ - void get(const double& posX, const double& posY, const double& posZ, - double& Bx, double& By, double& Bz) const override { - dd4hep::Position pos(posX,posY,posZ); - auto field = m_detector->field().magneticField(pos) * (dd4hep::kilogauss / dd4hep::tesla); - Bx = field.x(); - By = field.y(); - Bz = field.z(); - //return {field.x(), field.y(),field.z()}; - } - }; - -public: - DataHandle m_inputHitCollection{"inputHitCollection", Gaudi::DataHandle::Reader, this}; - DataHandle m_initialTrackParameters{"initialTrackParameters", Gaudi::DataHandle::Reader, this}; - DataHandle m_inputProtoTracks{"inputProtoTracks", Gaudi::DataHandle::Reader, this}; - DataHandle m_foundTracks{"trackParameters", Gaudi::DataHandle::Reader, this}; - DataHandle m_outputTrajectories{"outputTrajectories", Gaudi::DataHandle::Writer, this}; - - SmartIF m_geoSvc; - // Acts::GeometryContext m_geoctx; - // Acts::CalibrationContext m_calibctx; - // Acts::MagneticFieldContext m_fieldctx; - - std::map> m_detPlaneMap; - std::map m_surfaceMap; - - GenFitTrackFitter(const std::string& name, ISvcLocator* svcLoc); - - StatusCode initialize() override; - StatusCode execute() override; - }; - - -} // namespace Jug::Reco - -#endif diff --git a/JugTrack/src/components/HoughTransformProtoTracks.cpp b/JugTrack/src/components/HoughTransformProtoTracks.cpp index 1816b73..04eac00 100644 --- a/JugTrack/src/components/HoughTransformProtoTracks.cpp +++ b/JugTrack/src/components/HoughTransformProtoTracks.cpp @@ -16,7 +16,7 @@ #include "Math/Vector3D.h" -#include "eicd/TrackerHitCollection.h" +#include "edm4eic/TrackerHitCollection.h" namespace Jug::Reco { @@ -26,7 +26,7 @@ namespace Jug::Reco { */ class HoughTransformProtoTracks : public GaudiAlgorithm { private: - DataHandle m_inputTrackerHits{"inputTrackerHits", Gaudi::DataHandle::Reader, this}; + DataHandle m_inputTrackerHits{"inputTrackerHits", Gaudi::DataHandle::Reader, this}; DataHandle m_outputProtoTracks{"outputProtoTracks", Gaudi::DataHandle::Writer, this}; public: @@ -43,7 +43,7 @@ class HoughTransformProtoTracks : public GaudiAlgorithm { StatusCode execute() override { // input collection - //const eicd::TrackerHitCollection* hits = m_inputTrackerHits.get(); + //const edm4eic::TrackerHitCollection* hits = m_inputTrackerHits.get(); // Create output collections //auto proto_tracks = m_outputProtoTracks.createAndPut(); diff --git a/JugTrack/src/components/ParticlesFromTrackFit.cpp b/JugTrack/src/components/ParticlesFromTrackFit.cpp index 46af359..4312b65 100644 --- a/JugTrack/src/components/ParticlesFromTrackFit.cpp +++ b/JugTrack/src/components/ParticlesFromTrackFit.cpp @@ -22,16 +22,16 @@ #include "Acts/EventData/MultiTrajectoryHelpers.hpp" // Event Model related classes -#include "eicd/ReconstructedParticleCollection.h" -#include "eicd/TrackerHitCollection.h" -#include "eicd/TrackParametersCollection.h" +#include "edm4eic/ReconstructedParticleCollection.h" +#include "edm4eic/TrackerHitCollection.h" +#include "edm4eic/TrackParametersCollection.h" #include "JugTrack/IndexSourceLink.hpp" #include "JugTrack/Track.hpp" #include "JugTrack/Trajectories.hpp" #include "Acts/Utilities/Helpers.hpp" -#include "eicd/vector_utils.h" +#include "edm4eic/vector_utils.h" #include @@ -44,8 +44,8 @@ namespace Jug::Reco { class ParticlesFromTrackFit : public GaudiAlgorithm { private: DataHandle m_inputTrajectories{"inputTrajectories", Gaudi::DataHandle::Reader, this}; - DataHandle m_outputParticles{"outputParticles", Gaudi::DataHandle::Writer, this}; - DataHandle m_outputTrackParameters{"outputTrackParameters", Gaudi::DataHandle::Writer, this}; + DataHandle m_outputParticles{"outputParticles", Gaudi::DataHandle::Writer, this}; + DataHandle m_outputTrackParameters{"outputTrackParameters", Gaudi::DataHandle::Writer, this}; public: // ill-formed: using GaudiAlgorithm::GaudiAlgorithm; @@ -116,31 +116,31 @@ namespace Jug::Reco { debug() << " chi2 = " << trajState.chi2Sum << endmsg; } - const eicd::Vector2f loc { - parameter[Acts::eBoundLoc0], - parameter[Acts::eBoundLoc1] + const decltype(edm4eic::TrackParametersData::loc) loc { + static_cast(parameter[Acts::eBoundLoc0]), + static_cast(parameter[Acts::eBoundLoc1]) }; - const eicd::Cov3f covMomentum { + const decltype(edm4eic::TrackParametersData::momentumError) momentumError { static_cast(covariance(Acts::eBoundTheta, Acts::eBoundTheta)), static_cast(covariance(Acts::eBoundPhi, Acts::eBoundPhi)), static_cast(covariance(Acts::eBoundQOverP, Acts::eBoundQOverP)), static_cast(covariance(Acts::eBoundTheta, Acts::eBoundPhi)), static_cast(covariance(Acts::eBoundTheta, Acts::eBoundQOverP)), static_cast(covariance(Acts::eBoundPhi, Acts::eBoundQOverP))}; - const eicd::Cov2f covPos { + const decltype(edm4eic::TrackParametersData::locError) locError { static_cast(covariance(Acts::eBoundLoc0, Acts::eBoundLoc0)), static_cast(covariance(Acts::eBoundLoc1, Acts::eBoundLoc1)), static_cast(covariance(Acts::eBoundLoc0, Acts::eBoundLoc1))}; const float timeError{sqrt(static_cast(covariance(Acts::eBoundTime, Acts::eBoundTime)))}; - eicd::TrackParameters pars{ + edm4eic::TrackParameters pars{ 0, // type: track head --> 0 loc, - covPos, + locError, static_cast(parameter[Acts::eBoundTheta]), static_cast(parameter[Acts::eBoundPhi]), static_cast(parameter[Acts::eBoundQOverP]), - covMomentum, + momentumError, static_cast(parameter[Acts::eBoundTime]), timeError, static_cast(boundParam.charge())}; @@ -173,7 +173,7 @@ namespace Jug::Reco { auto rec_part = rec_parts->create(); rec_part.setMomentum( - eicd::sphericalToVector( + edm4eic::sphericalToVector( 1.0 / std::abs(params[Acts::eBoundQOverP]), params[Acts::eBoundTheta], params[Acts::eBoundPhi]) diff --git a/JugTrack/src/components/ProtoTrackMatching.cpp b/JugTrack/src/components/ProtoTrackMatching.cpp index 5065a54..5d4e3d9 100644 --- a/JugTrack/src/components/ProtoTrackMatching.cpp +++ b/JugTrack/src/components/ProtoTrackMatching.cpp @@ -16,7 +16,7 @@ #include "Math/Vector3D.h" -#include "eicd/TrackerHitCollection.h" +#include "edm4eic/TrackerHitCollection.h" namespace Jug::Reco { @@ -26,7 +26,7 @@ namespace Jug::Reco { */ class ProtoTrackMatching : public GaudiAlgorithm { private: - DataHandle m_inputTrackerHits{"inputTrackerHits", Gaudi::DataHandle::Reader, this}; + DataHandle m_inputTrackerHits{"inputTrackerHits", Gaudi::DataHandle::Reader, this}; DataHandle m_initialTrackParameters{"initialTrackParameters", Gaudi::DataHandle::Reader, this}; DataHandle m_inputProtoTracks{"inputProtoTracks", Gaudi::DataHandle::Reader, this}; DataHandle m_outputProtoTracks{"matchedProtoTracks", Gaudi::DataHandle::Writer, this}; diff --git a/JugTrack/src/components/SingleTrackSourceLinker.cpp b/JugTrack/src/components/SingleTrackSourceLinker.cpp index 27306f2..24fab9c 100644 --- a/JugTrack/src/components/SingleTrackSourceLinker.cpp +++ b/JugTrack/src/components/SingleTrackSourceLinker.cpp @@ -30,7 +30,7 @@ #include "JugTrack/Measurement.hpp" #include "JugTrack/ProtoTrack.hpp" -#include "eicd/TrackerHitCollection.h" +#include "edm4eic/TrackerHitCollection.h" namespace Jug::Reco { @@ -42,7 +42,7 @@ namespace Jug::Reco { */ class SingleTrackSourceLinker : public GaudiAlgorithm { private: - DataHandle m_inputHitCollection{"inputHitCollection", Gaudi::DataHandle::Reader, this}; + DataHandle m_inputHitCollection{"inputHitCollection", Gaudi::DataHandle::Reader, this}; DataHandle> m_sourceLinkStorage{"sourceLinkStorage", Gaudi::DataHandle::Writer, this}; DataHandle m_outputSourceLinks{"outputSourceLinks", Gaudi::DataHandle::Writer, this}; DataHandle m_outputMeasurements{"outputMeasurements", Gaudi::DataHandle::Writer, this}; @@ -74,7 +74,7 @@ class SingleTrackSourceLinker : public GaudiAlgorithm { StatusCode execute() override { // input collection - const eicd::TrackerHitCollection* hits = m_inputHitCollection.get(); + const edm4eic::TrackerHitCollection* hits = m_inputHitCollection.get(); // Create output collections auto* linkStorage = m_sourceLinkStorage.createAndPut(); auto* sourceLinks = m_outputSourceLinks.createAndPut(); diff --git a/JugTrack/src/components/TrackFindingAlgorithm.cpp b/JugTrack/src/components/TrackFindingAlgorithm.cpp index acdb678..1961d50 100644 --- a/JugTrack/src/components/TrackFindingAlgorithm.cpp +++ b/JugTrack/src/components/TrackFindingAlgorithm.cpp @@ -40,7 +40,7 @@ #include "JugTrack/Track.hpp" -#include "eicd/TrackerHitCollection.h" +#include "edm4eic/TrackerHitCollection.h" #include #include diff --git a/JugTrack/src/components/TrackFindingAlgorithm.h b/JugTrack/src/components/TrackFindingAlgorithm.h index dca6bce..b5eeca4 100644 --- a/JugTrack/src/components/TrackFindingAlgorithm.h +++ b/JugTrack/src/components/TrackFindingAlgorithm.h @@ -23,7 +23,7 @@ #include "JugTrack/Track.hpp" #include "JugTrack/Trajectories.hpp" -#include "eicd/TrackerHitCollection.h" +#include "edm4eic/TrackerHitCollection.h" #include "Acts/Definitions/Common.hpp" #include "Acts/Geometry/TrackingGeometry.hpp" diff --git a/JugTrack/src/components/TrackFittingAlgorithm.cpp b/JugTrack/src/components/TrackFittingAlgorithm.cpp index 6fa1344..644669f 100644 --- a/JugTrack/src/components/TrackFittingAlgorithm.cpp +++ b/JugTrack/src/components/TrackFittingAlgorithm.cpp @@ -36,7 +36,7 @@ #include "JugTrack/Track.hpp" #include "JugTrack/Measurement.hpp" -#include "eicd/TrackerHitCollection.h" +#include "edm4eic/TrackerHitCollection.h" #include #include diff --git a/JugTrack/src/components/TrackFittingAlgorithm.h b/JugTrack/src/components/TrackFittingAlgorithm.h index 6f193da..297b63d 100644 --- a/JugTrack/src/components/TrackFittingAlgorithm.h +++ b/JugTrack/src/components/TrackFittingAlgorithm.h @@ -24,7 +24,7 @@ #include "JugTrack/Trajectories.hpp" #include "JugTrack/ProtoTrack.hpp" -#include "eicd/TrackerHitCollection.h" +#include "edm4eic/TrackerHitCollection.h" #include "Acts/Definitions/Common.hpp" #include "Acts/Geometry/TrackingGeometry.hpp" diff --git a/JugTrack/src/components/TrackParamACTSSeeding.cpp b/JugTrack/src/components/TrackParamACTSSeeding.cpp index 6a664e6..8f302a9 100644 --- a/JugTrack/src/components/TrackParamACTSSeeding.cpp +++ b/JugTrack/src/components/TrackParamACTSSeeding.cpp @@ -31,7 +31,7 @@ #include "JugTrack/Measurement.hpp" #include "JugTrack/Track.hpp" -#include "eicd/TrackerHitCollection.h" +#include "edm4eic/TrackerHitCollection.h" #include "Math/Vector3D.h" @@ -63,7 +63,7 @@ namespace Jug::Reco { DataHandle m_inputMeasurements { "inputMeasurements", Gaudi::DataHandle::Reader, this}; - DataHandle + DataHandle m_inputHitCollection { "inputHitCollection", Gaudi::DataHandle::Reader, this }; DataHandle @@ -86,12 +86,12 @@ namespace Jug::Reco { /// Space point representation of eic::TrackerHitData suitable /// for ACTS track seeding. - class SpacePoint : eicd::TrackerHitData { + class SpacePoint : edm4eic::TrackerHitData { public: int32_t _measurementIndex; // Constructor to circumvent the fact that eic::TrackerHit // and associated classes are all non-polymorphic - SpacePoint(const eicd::TrackerHit h, + SpacePoint(const edm4eic::TrackerHit h, const int32_t measurementIndex) : _measurementIndex(measurementIndex) { @@ -212,7 +212,7 @@ namespace Jug::Reco { void findSeed(SeedContainer &seeds, - const eicd::TrackerHitCollection *hits, + const edm4eic::TrackerHitCollection *hits, const IndexSourceLinkContainer *sourceLinks, const MeasurementContainer *measurements, Acts::Seedfinder::State &state); @@ -291,7 +291,7 @@ namespace Jug::Reco { void TrackParamACTSSeeding:: findSeed(SeedContainer &seeds, - const eicd::TrackerHitCollection *hits, + const edm4eic::TrackerHitCollection *hits, const IndexSourceLinkContainer *sourceLinks, const MeasurementContainer *measurements, Acts::Seedfinder::State &state) @@ -330,11 +330,11 @@ namespace Jug::Reco { #ifdef USE_LOCAL_COORD spacePoint.push_back( SpacePoint( - eicd::TrackerHit( + edm4eic::TrackerHit( static_cast(spacePoint.size()), - eicd::Vector3f(v[0], v[1], v[2]), - eicd::CovDiag3f(25.0e-6 / 3.0, - 25.0e-6 / 3.0, 0.0), + {v[0], v[1], v[2]}, + {25.0e-6 / 3.0, + 25.0e-6 / 3.0, 0.0}, 0.0, 10.0, 0.05, 0.0), static_cast(spacePoint.size()))); spacePointPtrs.push_back(&spacePoint.back()); @@ -443,7 +443,7 @@ namespace Jug::Reco { StatusCode TrackParamACTSSeeding::execute() { - const eicd::TrackerHitCollection *hits = + const edm4eic::TrackerHitCollection *hits = m_inputHitCollection.get(); const IndexSourceLinkContainer *sourceLinks = m_inputSourceLinks.get(); diff --git a/JugTrack/src/components/TrackParamClusterInit.cpp b/JugTrack/src/components/TrackParamClusterInit.cpp index 538bcc8..e554a47 100644 --- a/JugTrack/src/components/TrackParamClusterInit.cpp +++ b/JugTrack/src/components/TrackParamClusterInit.cpp @@ -16,9 +16,9 @@ #include "JugBase/IGeoSvc.h" #include "JugTrack/Track.hpp" -#include "eicd/ClusterCollection.h" -#include "eicd/TrackerHitCollection.h" -#include "eicd/vector_utils.h" +#include "edm4eic/ClusterCollection.h" +#include "edm4eic/TrackerHitCollection.h" +#include "edm4eic/vector_utils.h" #include "Acts/Surfaces/PerigeeSurface.hpp" @@ -44,7 +44,7 @@ namespace Jug::Reco { */ class TrackParamClusterInit : public GaudiAlgorithm { private: - using Clusters = eicd::ClusterCollection; + using Clusters = edm4eic::ClusterCollection; DataHandle m_inputClusters{"inputClusters", Gaudi::DataHandle::Reader, this}; DataHandle m_outputInitialTrackParameters{"outputInitialTrackParameters", @@ -84,14 +84,14 @@ class TrackParamClusterInit : public GaudiAlgorithm { if (p < 0.1 * GeV) { continue; } - double len = eicd::magnitude(c.getPosition()); + double len = edm4eic::magnitude(c.getPosition()); auto momentum = c.getPosition() * p / len; Acts::BoundVector params; params(Acts::eBoundLoc0) = 0.0 * mm; params(Acts::eBoundLoc1) = 0.0 * mm; - params(Acts::eBoundPhi) = eicd::angleAzimuthal(momentum); - params(Acts::eBoundTheta) = eicd::anglePolar(momentum); + params(Acts::eBoundPhi) = edm4eic::angleAzimuthal(momentum); + params(Acts::eBoundTheta) = edm4eic::anglePolar(momentum); params(Acts::eBoundQOverP) = 1 / p; params(Acts::eBoundTime) = 0 * ns; @@ -107,8 +107,8 @@ class TrackParamClusterInit : public GaudiAlgorithm { Acts::BoundVector params2; params2(Acts::eBoundLoc0) = 0.0 * mm; params2(Acts::eBoundLoc1) = 0.0 * mm; - params2(Acts::eBoundPhi) = eicd::angleAzimuthal(momentum); - params2(Acts::eBoundTheta) = eicd::anglePolar(momentum); + params2(Acts::eBoundPhi) = edm4eic::angleAzimuthal(momentum); + params2(Acts::eBoundTheta) = edm4eic::anglePolar(momentum); params2(Acts::eBoundQOverP) = -1 / p; params2(Acts::eBoundTime) = 0 * ns; init_trk_params->push_back({pSurface, params2, -1}); diff --git a/JugTrack/src/components/TrackParamImagingClusterInit.cpp b/JugTrack/src/components/TrackParamImagingClusterInit.cpp index 7260951..5406384 100644 --- a/JugTrack/src/components/TrackParamImagingClusterInit.cpp +++ b/JugTrack/src/components/TrackParamImagingClusterInit.cpp @@ -16,9 +16,9 @@ #include "Acts/Definitions/Units.hpp" #include "Acts/Definitions/Common.hpp" -#include "eicd/TrackerHitCollection.h" -#include "eicd/ClusterCollection.h" -#include "eicd/vector_utils.h" +#include "edm4eic/TrackerHitCollection.h" +#include "edm4eic/ClusterCollection.h" +#include "edm4eic/vector_utils.h" #include "Acts/Surfaces/PerigeeSurface.hpp" @@ -44,7 +44,7 @@ namespace Jug::Reco { */ class TrackParamImagingClusterInit : public GaudiAlgorithm { private: - using ImagingClusters = eicd::ClusterCollection; + using ImagingClusters = edm4eic::ClusterCollection; DataHandle m_inputClusters{"inputClusters", Gaudi::DataHandle::Reader, this}; DataHandle m_outputInitialTrackParameters{"outputInitialTrackParameters", @@ -86,8 +86,8 @@ namespace Jug::Reco { if( p < 0.1*GeV) { continue; } - const double theta = eicd::anglePolar(c.getPosition()); - const double phi = eicd::angleAzimuthal(c.getPosition()); + const double theta = edm4eic::anglePolar(c.getPosition()); + const double phi = edm4eic::angleAzimuthal(c.getPosition()); Acts::BoundVector params; params(Acts::eBoundLoc0) = 0.0 * mm ; diff --git a/JugTrack/src/components/TrackParamTruthInit.cpp b/JugTrack/src/components/TrackParamTruthInit.cpp index 07b58b0..a07bd65 100644 --- a/JugTrack/src/components/TrackParamTruthInit.cpp +++ b/JugTrack/src/components/TrackParamTruthInit.cpp @@ -18,7 +18,7 @@ #include "Acts/Definitions/Units.hpp" #include "Acts/Definitions/Common.hpp" -#include "eicd/TrackerHitCollection.h" +#include "edm4eic/TrackerHitCollection.h" #include "edm4hep/MCParticleCollection.h" #include "Math/Vector3D.h" #include "Acts/Surfaces/PerigeeSurface.hpp" diff --git a/JugTrack/src/components/TrackParamVertexClusterInit.cpp b/JugTrack/src/components/TrackParamVertexClusterInit.cpp index 1e12819..362ce3b 100644 --- a/JugTrack/src/components/TrackParamVertexClusterInit.cpp +++ b/JugTrack/src/components/TrackParamVertexClusterInit.cpp @@ -19,9 +19,9 @@ #include "JugTrack/Track.hpp" #include "Acts/Surfaces/PerigeeSurface.hpp" -#include "eicd/ClusterCollection.h" -#include "eicd/TrackerHitCollection.h" -#include "eicd/vector_utils.h" +#include "edm4eic/ClusterCollection.h" +#include "edm4eic/TrackerHitCollection.h" +#include "edm4eic/vector_utils.h" using namespace Gaudi::Units; @@ -37,8 +37,8 @@ namespace Jug::Reco { */ class TrackParamVertexClusterInit : public GaudiAlgorithm { private: - using Clusters = eicd::ClusterCollection; - using VertexHits = eicd::TrackerHitCollection; + using Clusters = edm4eic::ClusterCollection; + using VertexHits = edm4eic::TrackerHitCollection; DataHandle m_inputVertexHits{"inputVertexHits", Gaudi::DataHandle::Reader, this}; DataHandle m_inputClusters{"inputClusters", Gaudi::DataHandle::Reader, this}; @@ -101,8 +101,8 @@ class TrackParamVertexClusterInit : public GaudiAlgorithm { Acts::BoundVector params; params(Acts::eBoundLoc0) = 0.0 * mm; params(Acts::eBoundLoc1) = 0.0 * mm; - params(Acts::eBoundPhi) = eicd::angleAzimuthal(momentum); - params(Acts::eBoundTheta) = eicd::anglePolar(momentum); + params(Acts::eBoundPhi) = edm4eic::angleAzimuthal(momentum); + params(Acts::eBoundTheta) = edm4eic::anglePolar(momentum); params(Acts::eBoundQOverP) = 1 / p_cluster; params(Acts::eBoundTime) = 0 * ns; @@ -114,8 +114,8 @@ class TrackParamVertexClusterInit : public GaudiAlgorithm { Acts::BoundVector params2; params2(Acts::eBoundLoc0) = 0.0 * mm; params2(Acts::eBoundLoc1) = 0.0 * mm; - params2(Acts::eBoundPhi) = eicd::angleAzimuthal(momentum); - params2(Acts::eBoundTheta) = eicd::anglePolar(momentum); + params2(Acts::eBoundPhi) = edm4eic::angleAzimuthal(momentum); + params2(Acts::eBoundTheta) = edm4eic::anglePolar(momentum); params2(Acts::eBoundQOverP) = -1 / p_cluster; params2(Acts::eBoundTime) = 0 * ns; init_trk_params->push_back({pSurface, params2, -1}); diff --git a/JugTrack/src/components/TrackProjector.cpp b/JugTrack/src/components/TrackProjector.cpp index 0a55fbf..620e1d1 100644 --- a/JugTrack/src/components/TrackProjector.cpp +++ b/JugTrack/src/components/TrackProjector.cpp @@ -22,10 +22,10 @@ #include "Acts/EventData/MultiTrajectoryHelpers.hpp" // Event Model related classes -#include "eicd/TrackerHitCollection.h" -#include "eicd/TrackParametersCollection.h" -#include "eicd/TrajectoryCollection.h" -#include "eicd/TrackSegmentCollection.h" +#include "edm4eic/TrackerHitCollection.h" +#include "edm4eic/TrackParametersCollection.h" +#include "edm4eic/TrajectoryCollection.h" +#include "edm4eic/TrackSegmentCollection.h" #include "JugTrack/IndexSourceLink.hpp" #include "JugTrack/Track.hpp" #include "JugTrack/Trajectories.hpp" @@ -38,7 +38,7 @@ #include "Acts/Propagator/EigenStepper.hpp" #include "Acts/Surfaces/PerigeeSurface.hpp" -#include "eicd/vector_utils.h" +#include "edm4eic/vector_utils.h" #include @@ -51,7 +51,7 @@ namespace Jug::Reco { class TrackProjector : public GaudiAlgorithm { private: DataHandle m_inputTrajectories{"inputTrajectories", Gaudi::DataHandle::Reader, this}; - DataHandle m_outputTrackSegments{"outputTrackSegments", Gaudi::DataHandle::Writer, this}; + DataHandle m_outputTrackSegments{"outputTrackSegments", Gaudi::DataHandle::Writer, this}; Gaudi::Property m_firstInVolumeID{this, "firstInVolumeID", 0}; Gaudi::Property m_firstInVolumeName{this, "firstInVolumeName", ""}; @@ -114,12 +114,10 @@ namespace Jug::Reco { debug() << "n state in trajectory " << m_nStates << endmsg; } - eicd::MutableTrackSegment track_segment; + edm4eic::MutableTrackSegment track_segment; // visit the track points mj.visitBackwards(trackTip, [&](auto&& trackstate) { - auto pathLength = trackstate.pathLength(); - // get volume info auto geoID = trackstate.referenceSurface().geometryId(); auto volume = geoID.volume(); @@ -139,29 +137,29 @@ namespace Jug::Reco { {0, 0, 0} ); // global position - const eicd::Vector3f position { - global.x(), - global.y(), - global.z() + const decltype(edm4eic::TrackPoint::position) position { + static_cast(global.x()), + static_cast(global.y()), + static_cast(global.z()) }; // local position - const eicd::Vector2f loc { - parameter[Acts::eBoundLoc0], - parameter[Acts::eBoundLoc1] + const decltype(edm4eic::TrackParametersData::loc) loc { + static_cast(parameter[Acts::eBoundLoc0]), + static_cast(parameter[Acts::eBoundLoc1]) }; - const eicd::Cov2f covLoc { + const decltype(edm4eic::TrackParametersData::locError) locError { static_cast(covariance(Acts::eBoundLoc0, Acts::eBoundLoc0)), static_cast(covariance(Acts::eBoundLoc1, Acts::eBoundLoc1)), static_cast(covariance(Acts::eBoundLoc0, Acts::eBoundLoc1)) }; - const eicd::Cov3f positionError{0, 0, 0}; - const eicd::Vector3f momentum = eicd::sphericalToVector( - 1.0 / std::abs(parameter[Acts::eBoundQOverP]), - parameter[Acts::eBoundTheta], - parameter[Acts::eBoundPhi] + const decltype(edm4eic::TrackPoint::positionError) positionError{0, 0, 0}; + const decltype(edm4eic::TrackPoint::momentum) momentum = edm4eic::sphericalToVector( + static_cast(1.0 / std::abs(parameter[Acts::eBoundQOverP])), + static_cast(parameter[Acts::eBoundTheta]), + static_cast(parameter[Acts::eBoundPhi]) ); - const eicd::Cov3f covMomentum { + const decltype(edm4eic::TrackPoint::momentumError) momentumError { static_cast(covariance(Acts::eBoundTheta, Acts::eBoundTheta)), static_cast(covariance(Acts::eBoundPhi, Acts::eBoundPhi)), static_cast(covariance(Acts::eBoundQOverP, Acts::eBoundQOverP)), @@ -173,11 +171,12 @@ namespace Jug::Reco { const float timeError{sqrt(static_cast(covariance(Acts::eBoundTime, Acts::eBoundTime)))}; const float theta(parameter[Acts::eBoundTheta]); const float phi(parameter[Acts::eBoundPhi]); - const eicd::Cov2f directionError { + const decltype(edm4eic::TrackPoint::directionError) directionError { static_cast(covariance(Acts::eBoundTheta, Acts::eBoundTheta)), static_cast(covariance(Acts::eBoundPhi, Acts::eBoundPhi)), static_cast(covariance(Acts::eBoundTheta, Acts::eBoundPhi)) }; + const float pathLength = static_cast(trackstate.pathLength()); const float pathLengthError = 0; // Store track point @@ -185,13 +184,13 @@ namespace Jug::Reco { position, positionError, momentum, - covMomentum, + momentumError, time, timeError, theta, phi, directionError, - static_cast(pathLength), + pathLength, pathLengthError }); diff --git a/JugTrack/src/components/TrackerSourceLinker.cpp b/JugTrack/src/components/TrackerSourceLinker.cpp index 4fd9cd8..3833d38 100644 --- a/JugTrack/src/components/TrackerSourceLinker.cpp +++ b/JugTrack/src/components/TrackerSourceLinker.cpp @@ -30,7 +30,7 @@ #include "JugTrack/IndexSourceLink.hpp" #include "JugTrack/Measurement.hpp" -#include "eicd/TrackerHitCollection.h" +#include "edm4eic/TrackerHitCollection.h" namespace Jug::Reco { @@ -44,7 +44,7 @@ namespace Jug::Reco { */ class TrackerSourceLinker : public GaudiAlgorithm { private: - DataHandle m_inputHitCollection{"inputHitCollection", Gaudi::DataHandle::Reader, this}; + DataHandle m_inputHitCollection{"inputHitCollection", Gaudi::DataHandle::Reader, this}; DataHandle> m_sourceLinkStorage{"sourceLinkStorage", Gaudi::DataHandle::Writer, this}; DataHandle m_outputSourceLinks{"outputSourceLinks", Gaudi::DataHandle::Writer, this}; DataHandle m_outputMeasurements{"outputMeasurements", Gaudi::DataHandle::Writer, this}; @@ -78,7 +78,7 @@ class TrackerSourceLinker : public GaudiAlgorithm { constexpr double mm_conv = mm_acts / dd4hep::mm; // = 1/0.1 // input collection - const eicd::TrackerHitCollection* hits = m_inputHitCollection.get(); + const edm4eic::TrackerHitCollection* hits = m_inputHitCollection.get(); // Create output collections auto* linkStorage = m_sourceLinkStorage.createAndPut(); auto* sourceLinks = m_outputSourceLinks.createAndPut(); diff --git a/JugTrack/src/components/TruthTrackSeeding.cpp b/JugTrack/src/components/TruthTrackSeeding.cpp index d650551..4a68a3a 100644 --- a/JugTrack/src/components/TruthTrackSeeding.cpp +++ b/JugTrack/src/components/TruthTrackSeeding.cpp @@ -20,7 +20,7 @@ //#include "Acts/EventData/Charge.hpp" #include "edm4hep/MCParticleCollection.h" -#include "eicd/TrackParametersCollection.h" +#include "edm4eic/TrackParametersCollection.h" #include "Math/Vector3D.h" @@ -28,7 +28,7 @@ namespace Jug::Reco { /** Track seeding using MC truth. * - * \note "Seeding" algorithms are required to output a eicd::TrackParametersCollection, as opposed to the legacy "init" + * \note "Seeding" algorithms are required to output a edm4eic::TrackParametersCollection, as opposed to the legacy "init" * algorithms, such as TrackParamTruthInit. * * \ingroup tracking @@ -37,7 +37,7 @@ namespace Jug::Reco { private: DataHandle m_inputMCParticles{"inputMCParticles", Gaudi::DataHandle::Reader, this}; - DataHandle m_outputTrackParameters{"outputTrackParameters", + DataHandle m_outputTrackParameters{"outputTrackParameters", Gaudi::DataHandle::Writer, this}; SmartIF m_pidSvc; @@ -93,7 +93,7 @@ namespace Jug::Reco { const auto q_over_p = charge / p; - eicd::TrackParameters params{-1, // type --> seed (-1) + edm4eic::TrackParameters params{-1, // type --> seed (-1) {0.0F, 0.0F}, // location on surface {0.1, 0.1, 0.1}, // Covariance on location theta, // theta (rad) diff --git a/doc/page.dox b/doc/page.dox index d0a904d..4a11fe2 100644 --- a/doc/page.dox +++ b/doc/page.dox @@ -22,7 +22,7 @@ * \brief Global geometry service * * \defgroup datasvc Data Services - * \brief Data services for using eicd/podio + * \brief Data services for using edm4eic/podio * */ diff --git a/external/algorithms/.clang-format b/external/algorithms/.clang-format new file mode 100644 index 0000000..5c443be --- /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 0000000..bb85179 --- /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/LICENSE b/external/algorithms/LICENSE new file mode 100644 index 0000000..c8fb052 --- /dev/null +++ b/external/algorithms/LICENSE @@ -0,0 +1,166 @@ + + GNU LESSER GENERAL PUBLIC LICENSE + Version 3, 29 June 2007 + + Copyright (C) 2007 Free Software Foundation, Inc. + Everyone is permitted to copy and distribute verbatim copies + of this license document, but changing it is not allowed. + + + This version of the GNU Lesser General Public License incorporates +the terms and conditions of version 3 of the GNU General Public +License, supplemented by the additional permissions listed below. + + 0. Additional Definitions. + + As used herein, "this License" refers to version 3 of the GNU Lesser +General Public License, and the "GNU GPL" refers to version 3 of the GNU +General Public License. + + "The Library" refers to a covered work governed by this License, +other than an Application or a Combined Work as defined below. + + An "Application" is any work that makes use of an interface provided +by the Library, but which is not otherwise based on the Library. +Defining a subclass of a class defined by the Library is deemed a mode +of using an interface provided by the Library. + + A "Combined Work" is a work produced by combining or linking an +Application with the Library. The particular version of the Library +with which the Combined Work was made is also called the "Linked +Version". + + The "Minimal Corresponding Source" for a Combined Work means the +Corresponding Source for the Combined Work, excluding any source code +for portions of the Combined Work that, considered in isolation, are +based on the Application, and not on the Linked Version. + + The "Corresponding Application Code" for a Combined Work means the +object code and/or source code for the Application, including any data +and utility programs needed for reproducing the Combined Work from the +Application, but excluding the System Libraries of the Combined Work. + + 1. Exception to Section 3 of the GNU GPL. + + You may convey a covered work under sections 3 and 4 of this License +without being bound by section 3 of the GNU GPL. + + 2. Conveying Modified Versions. + + If you modify a copy of the Library, and, in your modifications, a +facility refers to a function or data to be supplied by an Application +that uses the facility (other than as an argument passed when the +facility is invoked), then you may convey a copy of the modified +version: + + a) under this License, provided that you make a good faith effort to + ensure that, in the event an Application does not supply the + function or data, the facility still operates, and performs + whatever part of its purpose remains meaningful, or + + b) under the GNU GPL, with none of the additional permissions of + this License applicable to that copy. + + 3. Object Code Incorporating Material from Library Header Files. + + The object code form of an Application may incorporate material from +a header file that is part of the Library. You may convey such object +code under terms of your choice, provided that, if the incorporated +material is not limited to numerical parameters, data structure +layouts and accessors, or small macros, inline functions and templates +(ten or fewer lines in length), you do both of the following: + + a) Give prominent notice with each copy of the object code that the + Library is used in it and that the Library and its use are + covered by this License. + + b) Accompany the object code with a copy of the GNU GPL and this license + document. + + 4. Combined Works. + + You may convey a Combined Work under terms of your choice that, +taken together, effectively do not restrict modification of the +portions of the Library contained in the Combined Work and reverse +engineering for debugging such modifications, if you also do each of +the following: + + a) Give prominent notice with each copy of the Combined Work that + the Library is used in it and that the Library and its use are + covered by this License. + + b) Accompany the Combined Work with a copy of the GNU GPL and this license + document. + + c) For a Combined Work that displays copyright notices during + execution, include the copyright notice for the Library among + these notices, as well as a reference directing the user to the + copies of the GNU GPL and this license document. + + d) Do one of the following: + + 0) Convey the Minimal Corresponding Source under the terms of this + License, and the Corresponding Application Code in a form + suitable for, and under terms that permit, the user to + recombine or relink the Application with a modified version of + the Linked Version to produce a modified Combined Work, in the + manner specified by section 6 of the GNU GPL for conveying + Corresponding Source. + + 1) Use a suitable shared library mechanism for linking with the + Library. A suitable mechanism is one that (a) uses at run time + a copy of the Library already present on the user's computer + system, and (b) will operate properly with a modified version + of the Library that is interface-compatible with the Linked + Version. + + e) Provide Installation Information, but only if you would otherwise + be required to provide such information under section 6 of the + GNU GPL, and only to the extent that such information is + necessary to install and execute a modified version of the + Combined Work produced by recombining or relinking the + Application with a modified version of the Linked Version. (If + you use option 4d0, the Installation Information must accompany + the Minimal Corresponding Source and Corresponding Application + Code. If you use option 4d1, you must provide the Installation + Information in the manner specified by section 6 of the GNU GPL + for conveying Corresponding Source.) + + 5. Combined Libraries. + + You may place library facilities that are a work based on the +Library side by side in a single library together with other library +facilities that are not Applications and are not covered by this +License, and convey such a combined library under terms of your +choice, if you do both of the following: + + a) Accompany the combined library with a copy of the same work based + on the Library, uncombined with any other library facilities, + conveyed under the terms of this License. + + b) Give prominent notice with the combined library that part of it + is a work based on the Library, and explaining where to find the + accompanying uncombined form of the same work. + + 6. Revised Versions of the GNU Lesser General Public License. + + The Free Software Foundation may publish revised and/or new versions +of the GNU Lesser General Public License from time to time. Such new +versions will be similar in spirit to the present version, but may +differ in detail to address new problems or concerns. + + Each version is given a distinguishing version number. If the +Library as you received it specifies that a certain numbered version +of the GNU Lesser General Public License "or any later version" +applies to it, you have the option of following the terms and +conditions either of that published version or of any later version +published by the Free Software Foundation. If the Library as you +received it does not specify a version number of the GNU Lesser +General Public License, you may choose any version of the GNU Lesser +General Public License ever published by the Free Software Foundation. + + If the Library as you received it specifies that a proxy can decide +whether future versions of the GNU Lesser General Public License shall +apply, that proxy's public statement of acceptance of any version is +permanent authorization for you to choose that version for the +Library. diff --git a/external/algorithms/LICENSE.spdx b/external/algorithms/LICENSE.spdx new file mode 100644 index 0000000..c76ebdc --- /dev/null +++ b/external/algorithms/LICENSE.spdx @@ -0,0 +1,6 @@ +SPDXVersion: SPDX-2.1 +DataLicense: CC0-1.0 +PackageName: Juggler +PackageOriginator: Wouter Deconinck, Sylvester Joosten +PackageHomePage: https://github.com/eic/algorithms +PackageLicenseDeclared: LGPL-3.0-or-later diff --git a/external/algorithms/calorimetry/CMakeLists.txt b/external/algorithms/calorimetry/CMakeLists.txt new file mode 100644 index 0000000..b235229 --- /dev/null +++ b/external/algorithms/calorimetry/CMakeLists.txt @@ -0,0 +1,45 @@ +# SPDX-License-Identifier: LGPL-3.0-or-later +# Copyright (C) 2022 Wouter Deconinck, 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/ClusterRecoCoG.h b/external/algorithms/calorimetry/include/algorithms/calorimetry/ClusterRecoCoG.h new file mode 100644 index 0000000..b43090a --- /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 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/ClusterRecoCoG.cpp b/external/algorithms/calorimetry/src/ClusterRecoCoG.cpp new file mode 100644 index 0000000..9501871 --- /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 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 0000000..2b5d16d --- /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 0000000..2591084 --- /dev/null +++ b/external/algorithms/core/CMakeLists.txt @@ -0,0 +1,44 @@ +# SPDX-License-Identifier: LGPL-3.0-or-later +# Copyright (C) 2022 Wouter Deconinck, 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 0000000..3fce167 --- /dev/null +++ b/external/algorithms/core/include/algorithms/algorithm.h @@ -0,0 +1,95 @@ +// SPDX-License-Identifier: LGPL-3.0-or-later +// Copyright (C) 2022 Wouter Deconinck, Sylvester Joosten +// +// Algorithm base class, defined as a template with a tuple of Input<> and Output<> +// parameters. +// +// Known types are: +// - Normal data type T +// - Optional data type std::optional +// - Vector of normal data std::vector +// +// For input data, this then selects: +// - T --> gsl::not_null (NOT allowed to be null) +// - optional --> const T* (allowed to be null) +// - vector --> std::vector> (N arguments, NOT allowed to +// be null, but can be zero +// length) +// +// Same for output data, but replace `const T*` with `T*` (mutable) everywhere. +// +// The ::process() algorithm is then provided with a tuple of both the input and the +// output pointers according to this scheme. +// +// Finally, provides provides utility traits to determine if a type Input or Output are +// an Input or Output Type (is_input_v and is_output_v) +// +#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 0000000..4d96e30 --- /dev/null +++ b/external/algorithms/core/include/algorithms/detail/upcast.h @@ -0,0 +1,36 @@ +#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 <> struct upcast_type { using type = bool; }; + +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 0000000..ced872b --- /dev/null +++ b/external/algorithms/core/include/algorithms/error.h @@ -0,0 +1,28 @@ +// SPDX-License-Identifier: LGPL-3.0-or-later +// Copyright (C) 2022 Wouter Deconinck, Sylvester Joosten +// +// A simple exception base class that transports both an error message and error type + +#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 0000000..7617796 --- /dev/null +++ b/external/algorithms/core/include/algorithms/geo.h @@ -0,0 +1,46 @@ +// SPDX-License-Identifier: LGPL-3.0-or-later +// Copyright (C) 2022 Wouter Deconinck, Sylvester Joosten +// +// DD4hep Geometry service exposing a detector(), world(), and cellIDPositionConverter() +// Meant to be set by the calling framework, but can also load DD4hep itself +// when given a compact file as Property. +// +#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 0000000..26799f7 --- /dev/null +++ b/external/algorithms/core/include/algorithms/logger.h @@ -0,0 +1,243 @@ +// SPDX-License-Identifier: LGPL-3.0-or-later +// Copyright (C) 2022 Wouter Deconinck, Sylvester Joosten +// +// Logging service, which will use a callback to use the framework logging infrastructure. +// use ::action(void(LogLevel, std::string_view, std::string_view)) to register +// a logger. +// +// Also provides the LoggerMixin and LoggedService base classes +// +#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 0000000..e80875e --- /dev/null +++ b/external/algorithms/core/include/algorithms/name.h @@ -0,0 +1,21 @@ +// SPDX-License-Identifier: LGPL-3.0-or-later +// Copyright (C) 2022 Wouter Deconinck, Sylvester Joosten +// +// Defines NameMixin - simple base class to provide name functionality +// +#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 0000000..02f2716 --- /dev/null +++ b/external/algorithms/core/include/algorithms/property.h @@ -0,0 +1,176 @@ +// SPDX-License-Identifier: LGPL-3.0-or-later +// Copyright (C) 2022 Wouter Deconinck, Sylvester Joosten +// +// Defines the Configurable base class and related PropertyMixin +// These base classes provide access to the Property object, a self-registering +// configurable property that acts as a bare object T from a performance point-of-view +// +#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/service.h b/external/algorithms/core/include/algorithms/service.h new file mode 100644 index 0000000..a6a5604 --- /dev/null +++ b/external/algorithms/core/include/algorithms/service.h @@ -0,0 +1,80 @@ +// SPDX-License-Identifier: LGPL-3.0-or-later +// Copyright (C) 2022 Wouter Deconinck, Sylvester Joosten +// +// Basic Service infrastructure, implemented as lazy-evaluated singletons (thread-safe). +// +// Provides the special ServiceSvc (that provides access to all instantiated services as +// Configurable*), and the ServiceMixin (and related ALGORITHMS_DEFINE_SERVICE macro). +// +#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 0000000..c7fe8a0 --- /dev/null +++ b/external/algorithms/core/include/algorithms/type_traits.h @@ -0,0 +1,53 @@ +// SPDX-License-Identifier: LGPL-3.0-or-later +// Copyright (C) 2022 Wouter Deconinck, Sylvester Joosten +// +// Type traits used for argument deduction for input and output arguments. +// It allows to distinguish Vector and Optional arguments from regular arguments, +// and to select the appropriate underlying pointer type for each of the argument types. +// +#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 0000000..a264023 --- /dev/null +++ b/external/algorithms/core/src/dummy.cpp @@ -0,0 +1,91 @@ +// SPDX-License-Identifier: LGPL-3.0-or-later +// Copyright (C) 2022 Wouter Deconinck, Sylvester Joosten +// +// Dummy that instantiates some objects to trigger compile errors early on in case of +// bugs. This should be migrated to tests TODO. +#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 0000000..e910767 --- /dev/null +++ b/external/algorithms/core/src/geo.cpp @@ -0,0 +1,32 @@ +// SPDX-License-Identifier: LGPL-3.0-or-later +// Copyright (C) 2022 Wouter Deconinck, Sylvester Joosten +// +// Implementation of the geo service +// +#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 +