From 7673cb84c0c3ec013da249bbca7e6feccf3af240 Mon Sep 17 00:00:00 2001 From: Dmitry Romanov Date: Sun, 10 Nov 2024 12:11:02 -0500 Subject: [PATCH 1/8] PadGeometryHelper better quality code --- source/tdis/PadGeometryHelper.hpp | 68 ++++++++++++++++--------------- 1 file changed, 35 insertions(+), 33 deletions(-) diff --git a/source/tdis/PadGeometryHelper.hpp b/source/tdis/PadGeometryHelper.hpp index ef9fb03..8d24484 100644 --- a/source/tdis/PadGeometryHelper.hpp +++ b/source/tdis/PadGeometryHelper.hpp @@ -13,46 +13,48 @@ constexpr double total_radius = max_radius - min_radius; // Total radial width constexpr double ring_width = total_radius / num_rings; // Radial width of each ring constexpr double delta_theta = 2 * M_PI / num_pads_per_ring; // Angular width of each pad in radians -inline std::pair getPadCenter(int ring, int pad) { - /* - Compute the X and Y coordinates of the center of a pad given its ring and pad indices. +namespace tdis { + inline std::tuple getPadCenter(int ring, int pad) { + /* + Compute the X and Y coordinates of the center of a pad given its ring and pad indices. - Parameters: - - ring (int): The ring index (0 is the innermost ring). - - pad (int): The pad index (0 is at or closest to φ=0°, numbering is clockwise). + Parameters: + - ring (int): The ring index (0 is the innermost ring). + - pad (int): The pad index (0 is at or closest to φ=0°, numbering is clockwise). - Returns: - - std::pair: X and Y coordinates of the pad center in cm. - */ + Returns: + - std::pair: X and Y coordinates of the pad center in cm. + */ - // Validate inputs - if (ring < 0 || ring >= num_rings) { - throw std::invalid_argument("Ring index must be between 0 and " + std::to_string(num_rings - 1)); - } - if (pad < 0 || pad >= num_pads_per_ring) { - throw std::invalid_argument("Pad index must be between 0 and " + std::to_string(num_pads_per_ring - 1)); - } + // Validate inputs + if (ring < 0 || ring >= num_rings) { + throw std::invalid_argument("Ring index must be between 0 and " + std::to_string(num_rings - 1)); + } + if (pad < 0 || pad >= num_pads_per_ring) { + throw std::invalid_argument("Pad index must be between 0 and " + std::to_string(num_pads_per_ring - 1)); + } - // Compute radial center of the ring - double r_center = min_radius + ring_width * (ring + 0.5); + // Compute radial center of the ring + const double r_center = min_radius + ring_width * (ring + 0.5); - // Determine the angular offset for odd rings - double theta_offset = (ring % 2 == 0) ? 0.0 : delta_theta / 2.0; + // Determine the angular offset for odd rings + const double theta_offset = (ring % 2 == 0) ? 0.0 : delta_theta / 2.0; - // Compute the angular center of the pad (in radians) - double theta_clockwise = pad * delta_theta + theta_offset + delta_theta / 2.0; + // Compute the angular center of the pad (in radians) + const double theta_clockwise = pad * delta_theta + theta_offset + delta_theta / 2.0; - // Convert to counterclockwise angle for standard coordinate system - double theta_center = 2 * M_PI - theta_clockwise; + // Convert to counterclockwise angle for standard coordinate system + double theta_center = 2 * M_PI - theta_clockwise; - // Ensure theta_center is within [0, 2π] - if (theta_center < 0) { - theta_center += 2 * M_PI; - } + // Ensure theta_center is within [0, 2π] + if (theta_center < 0) { + theta_center += 2 * M_PI; + } - // Convert from polar to Cartesian coordinates - double x = r_center * std::cos(theta_center); - double y = r_center * std::sin(theta_center); + // Convert from polar to Cartesian coordinates + double x = r_center * std::cos(theta_center); + double y = r_center * std::sin(theta_center); - return { x, y }; -} + return { x, y }; + } +} \ No newline at end of file From 95895a1b9e0f755dc05175a6cd3fe61565341e7e Mon Sep 17 00:00:00 2001 From: Dmitry Romanov Date: Mon, 11 Nov 2024 16:34:15 -0500 Subject: [PATCH 2/8] Refactor data model, cleanup cmake, progress tracking --- source/tdis/CMakeLists.txt | 74 +++---------------- source/tdis/io/DigitizedDataEventSource.hpp | 8 +- source/tdis/io/PodioWriteProcessor.hpp | 42 ++++++++--- source/tdis/layout.yaml | 6 +- source/tdis/services/LogService.hpp | 2 +- source/tdis/tdis_main.cpp | 4 +- .../{ => tracking}/CKFTrackingFunction.cc | 0 source/tdis/tracking/CkfFitFactory.h | 36 +++++++++ 8 files changed, 88 insertions(+), 84 deletions(-) rename source/tdis/{ => tracking}/CKFTrackingFunction.cc (100%) create mode 100644 source/tdis/tracking/CkfFitFactory.h diff --git a/source/tdis/CMakeLists.txt b/source/tdis/CMakeLists.txt index d5c0a18..62765f0 100644 --- a/source/tdis/CMakeLists.txt +++ b/source/tdis/CMakeLists.txt @@ -13,6 +13,7 @@ set(CMAKE_CXX_EXTENSIONS OFF) cmake_policy(SET CMP0144 NEW) +# --------- FETCH DEPENDENCIES -------------- # Fetch CLI11 include(FetchContent) FetchContent_Declare( @@ -23,13 +24,11 @@ FetchContent_Declare( FetchContent_MakeAvailable(CLI11) - # --------- GENERATE PODIO -------------- find_package(podio REQUIRED) PODIO_GENERATE_DATAMODEL(podio_model "layout.yaml" DATAMODEL_HEADERS DATAMODEL_SOURCES OUTPUT_FOLDER "${CMAKE_CURRENT_SOURCE_DIR}/podio_model") PODIO_ADD_DATAMODEL_CORE_LIB(podio_model_lib "${DATAMODEL_HEADERS}" "${DATAMODEL_SOURCES}") PODIO_ADD_ROOT_IO_DICT(podio_model_dict podio_model_lib "${DATAMODEL_HEADERS}" "podio_model/src/selection.xml") -message(STATUS "PODIO MODEL DICT: ${podio_model_dict}") target_include_directories(podio_model_lib PRIVATE "podio_model") target_include_directories(podio_model_dict PRIVATE "podio_model") @@ -55,11 +54,20 @@ find_package(JANA REQUIRED) find_package(fmt REQUIRED) find_package(spdlog REQUIRED) find_package(Boost REQUIRED) - +find_package(ROOT COMPONENTS Core RIO Hist Graf Gpad Tree Postscript Matrix Physics MathCore) find_package(Acts REQUIRED COMPONENTS Core PluginTGeo PluginJson) target_include_directories(tdis PUBLIC ${CMAKE_CURRENT_LIST_DIR} "${CMAKE_CURRENT_LIST_DIR}/..") -target_link_libraries(tdis ${JANA_LIB} podio::podio podio::podioRootIO podio_model_lib podio_model_dict spdlog::spdlog fmt::fmt CLI11::CLI11 Boost::boost ActsCore ActsPluginTGeo ActsPluginJson) +target_link_libraries(tdis + ${JANA_LIB} + ROOT::RIO ROOT::Core + podio::podio podio::podioRootIO podio_model_lib podio_model_dict + spdlog::spdlog + fmt::fmt + CLI11::CLI11 + Boost::boost + ActsCore ActsPluginTGeo ActsPluginJson +) set_target_properties(tdis PROPERTIES INSTALL_RPATH_USE_LINK_PATH TRUE) target_include_directories(tdis PRIVATE "podio_model") @@ -103,61 +111,3 @@ if(WITH_TESTS) message(WARNING "Catch2 not found, unit tests will not be built.") endif() endif() - -# -# -## Function creates ${PLUGIN_NAME}_plugin and ${PLUGIN_NAME}_library targets -## Setting default includes, libraries and installation paths -#plugin_add(${PLUGIN_NAME}) -# -## The macro grabs sources as *.cc *.cpp *.c and headers as *.h *.hh *.hpp Then -## correctly sets sources for ${_name}_plugin and ${_name}_library targets Adds -## headers to the correct installation directory -#plugin_glob_all(${PLUGIN_NAME}) -# -## Find dependencies -# -#plugin_add_acts(${PLUGIN_NAME}) -#plugin_add_cern_root(${PLUGIN_NAME}) -#plugin_add_event_model(${PLUGIN_NAME}) -# -# -# -#plugin_sources(${PLUGIN_NAME} ${DATAMODEL_HEADERS} ${DATAMODEL_SOURCES}) - - -#set(PodioExample_SOURCES -# PodioExample.cc -# PodioExampleProcessor.cc -# PodioExampleSource.cc -# ExampleClusterFactory.cc -# ExampleMultifactory.cc -# ) -# -#foreach( _conf ${CMAKE_CONFIGURATION_TYPES} ) -# string(TOUPPER ${_conf} _conf ) -# set( CMAKE_RUNTIME_OUTPUT_DIRECTORY_${_conf} ${CMAKE_CURRENT_BINARY_DIR} ) -# set( CMAKE_LIBRARY_OUTPUT_DIRECTORY_${_conf} ${CMAKE_CURRENT_BINARY_DIR} ) -# set( CMAKE_ARCHIVE_OUTPUT_DIRECTORY_${_conf} ${CMAKE_CURRENT_BINARY_DIR} ) -#endforeach() -# -#PODIO_GENERATE_DATAMODEL(datamodel layout.yaml DATAMODEL_HEADERS DATAMODEL_SOURCES IO_BACKEND_HANDLERS ROOT) -# -#PODIO_ADD_DATAMODEL_CORE_LIB(PodioExampleDatamodel "${DATAMODEL_HEADERS}" "${DATAMODEL_SOURCES}") -# -#PODIO_ADD_ROOT_IO_DICT(PodioExampleDatamodelDict PodioExampleDatamodel "${DATAMODEL_HEADERS}" src/selection.xml) -# -# -# -#add_executable(PodioExample ${PodioExample_SOURCES}) -#target_include_directories(PodioExample PUBLIC .) -#target_link_libraries(PodioExample jana2 podio::podio PodioExampleDatamodel PodioExampleDatamodelDict podio::podioRootIO) -#set_target_properties(PodioExample PROPERTIES INSTALL_RPATH_USE_LINK_PATH TRUE) -# -#install(TARGETS PodioExample DESTINATION bin) -#install(TARGETS PodioExampleDatamodel DESTINATION lib) -#install(TARGETS PodioExampleDatamodelDict DESTINATION lib) - - - - diff --git a/source/tdis/io/DigitizedDataEventSource.hpp b/source/tdis/io/DigitizedDataEventSource.hpp index b0b0906..589ab40 100644 --- a/source/tdis/io/DigitizedDataEventSource.hpp +++ b/source/tdis/io/DigitizedDataEventSource.hpp @@ -55,7 +55,7 @@ #include "podio_model/EventInfoCollection.h" #include "services/LogService.hpp" -namespace tdis::readout { +namespace tdis::io { /** POD structure for readout hits **/ struct DigitizedReadoutHit { @@ -336,8 +336,8 @@ namespace tdis::readout { EventInfoCollection info; info.push_back(MutableEventInfo(0, 0, 0)); // event nr, timeslice nr, run nr event.InsertCollection(std::move(info), "EventInfo"); - event.InsertCollection(std::move(podioTracks), "McTracks"); - event.InsertCollection(std::move(podioHits), "McHits"); + event.InsertCollection(std::move(podioTracks), "DigitizedMtpcMcTrack"); + event.InsertCollection(std::move(podioHits), "DigitizedMtpcMcHit"); m_log->info("Event has been emitted at {}", m_event_line_index); return Result::Success; } @@ -352,7 +352,7 @@ namespace tdis::readout { // The template specialization needs to be in the global namespace (or at least not inside the tdis namespace) template <> -inline double JEventSourceGeneratorT::CheckOpenable(std::string resource_name) { +inline double JEventSourceGeneratorT::CheckOpenable(std::string resource_name) { // CheckOpenable() decides how confident we are that this EventSource can handle this resource. // 0.0 -> 'Cannot handle' // (0.0, 1.0] -> 'Cean handle, with this confidence level' diff --git a/source/tdis/io/PodioWriteProcessor.hpp b/source/tdis/io/PodioWriteProcessor.hpp index 90c1af2..c060335 100644 --- a/source/tdis/io/PodioWriteProcessor.hpp +++ b/source/tdis/io/PodioWriteProcessor.hpp @@ -64,7 +64,7 @@ #include "services/LogService.hpp" -//namespace tdis::readout { +namespace tdis::io { class PodioWriteProcessor : public JEventProcessor { public: @@ -74,8 +74,8 @@ class PodioWriteProcessor : public JEventProcessor { "EventInfo", // Truth record - "McTracks", - "McHits" + "DigitizedMtpcMcTrack", + "DigitizedMtpcMcHit" }; PodioWriteProcessor(JApplication * app); @@ -127,25 +127,36 @@ inline void PodioWriteProcessor::Init() { "be written (including ones from input file). Don't set this and use " "PODIO:OUTPUT_EXCLUDE_COLLECTIONS to write everything except a selection."); - m_output_collections = std::set(output_collections.begin(), output_collections.end()); + m_output_collections + = std::set(output_collections.begin(), output_collections.end()); m_app->SetDefaultParameter( - "podio:print_collections", - m_collections_to_print, - "Comma separated list of collection names to print to screen, e.g. for debugging." - ); + "podio:print_collections", m_collections_to_print, + "Comma separated list of collection names to print to screen, e.g. for debugging."); m_writer = std::make_unique(m_output_file); } - -void PodioWriteProcessor::Process(const std::shared_ptr& event) { +inline void PodioWriteProcessor::Process(const std::shared_ptr& event) { std::lock_guard lock(m_mutex); + m_log->info("PodioWriteProcessor::Process() All event collections:"); + auto event_collections = event->GetAllCollectionNames(); + for (const auto& coll_name : event_collections) { + try { + m_log->info(" {}", coll_name); + } catch (std::exception& e) { + // chomp + } + } + + // Trigger all collections once to fix the collection IDs + m_collections_to_write.clear(); for (const auto& coll_name : m_output_collections) { try { [[maybe_unused]] const auto* coll_ptr = event->GetCollectionBase(coll_name); + m_collections_to_write.push_back(coll_name); } catch (std::exception& e) { // chomp } @@ -243,11 +254,18 @@ void PodioWriteProcessor::Process(const std::shared_ptr& event) { } */ m_writer->writeFrame(*frame, "events", m_collections_to_write); + + auto [missing_names, all_names] = m_writer->checkConsistency(m_collections_to_write, ""); + m_log->info("PODIO checkConsistency missing_names: {}", missing_names.size()); + for (const auto& coll_name : missing_names) { + m_log->info(" {}", coll_name); + + } m_is_first_event = false; } -void PodioWriteProcessor::Finish() { +inline void PodioWriteProcessor::Finish() { m_writer->finish(); } -//} \ No newline at end of file +} // namespace tdis:io \ No newline at end of file diff --git a/source/tdis/layout.yaml b/source/tdis/layout.yaml index b2c3731..a63ca6f 100644 --- a/source/tdis/layout.yaml +++ b/source/tdis/layout.yaml @@ -252,7 +252,7 @@ components: - std::array transform // row-wise 4x4 affine transform [R T; 0 1] with 3x3 rotation matrix R and translation column 3-vector T datatypes : - EventInfo: + tdis::EventInfo: Description : "Event info" Author : "N. Brei" Members : @@ -260,14 +260,14 @@ datatypes : - int TimesliceNumber // timeslice number - int RunNumber // run number - TimesliceInfo: + tdis::TimesliceInfo: Description : "Timeslice info" Author : "N. Brei" Members : - int TimesliceNumber // timeslice number - int RunNumber // run number - DigitizedMtpcMcTrack: + tdis::DigitizedMtpcMcTrack: Description: "TDIS MTPC Digitized track" Author: "Dmitry Romanov" Members: diff --git a/source/tdis/services/LogService.hpp b/source/tdis/services/LogService.hpp index 54b8493..9f5035f 100644 --- a/source/tdis/services/LogService.hpp +++ b/source/tdis/services/LogService.hpp @@ -111,7 +111,7 @@ namespace tdis::services { // Set log level for this named logger allowing user to specify as config. parameter // e.g. EcalEndcapPRecHits:LogLevel std::string log_level_str = default_level ? LogLevelToString(default_level.value()) : m_log_level_str; - m_application->SetDefaultParameter(name + ":LogLevel", log_level_str, "log_level for " + name + ": trace, debug, info, warn, error, critical, off"); + m_application->SetDefaultParameter(name + ":log_level", log_level_str, "log_level for " + name + ": trace, debug, info, warn, error, critical, off"); logger->set_level(ParseLogLevel(log_level_str)); } return logger; diff --git a/source/tdis/tdis_main.cpp b/source/tdis/tdis_main.cpp index 3156091..6fc9c81 100644 --- a/source/tdis/tdis_main.cpp +++ b/source/tdis/tdis_main.cpp @@ -109,8 +109,8 @@ int main(int argc, char* argv[]) { JApplication app(parameterManager); - app.Add(new JEventSourceGeneratorT); - app.Add(new PodioWriteProcessor(&app)); + app.Add(new JEventSourceGeneratorT); + app.Add(new tdis::io::PodioWriteProcessor(&app)); app.ProvideService(std::make_shared(&app)); // app.Add(new JEventProcessorPodio); // app.Add(new JFactoryGeneratorT()); diff --git a/source/tdis/CKFTrackingFunction.cc b/source/tdis/tracking/CKFTrackingFunction.cc similarity index 100% rename from source/tdis/CKFTrackingFunction.cc rename to source/tdis/tracking/CKFTrackingFunction.cc diff --git a/source/tdis/tracking/CkfFitFactory.h b/source/tdis/tracking/CkfFitFactory.h new file mode 100644 index 0000000..f6c99a6 --- /dev/null +++ b/source/tdis/tracking/CkfFitFactory.h @@ -0,0 +1,36 @@ +#pragma once + +#include +#include + +namespace tdis::tracking { + + + + struct MyClusterFactory : public JOmniFactory { + + PodioInput m_protoclusters_in {this}; + PodioOutput m_clusters_out {this}; + + + void Configure() { + } + + void ChangeRun(int32_t /*run_nr*/) { + } + + void Execute(int32_t /*run_nr*/, uint64_t /*evt_nr*/) { + + auto cs = std::make_unique(); + + for (auto protocluster : *m_protoclusters_in()) { + auto cluster = MutableExampleCluster(protocluster.energy() + 1000); + cluster.addClusters(protocluster); + cs->push_back(cluster); + } + + m_clusters_out() = std::move(cs); + } + }; + +} \ No newline at end of file From 0c445411ba124934bd953bb524cff8f43cb85c3b Mon Sep 17 00:00:00 2001 From: Dmitry Romanov Date: Thu, 14 Nov 2024 22:05:23 -0500 Subject: [PATCH 3/8] Geometry and Data source updates --- docker/tdis-pre/Dockerfile | 4 +- source/tdis/CKFTracking.h | 8 +- source/tdis/CMakeLists.txt | 40 +++ source/tdis/io/DigitizedDataEventSource.hpp | 46 ++- source/tdis/layout.yaml | 3 + source/tdis/tdis_main.cpp | 11 +- source/tdis/tracking/ActsGeometryService.cc | 327 +++++++++++++++++- source/tdis/tracking/ActsGeometryService.h | 64 ++-- .../tdis/tracking/BuildTelescopeDetector.cpp | 131 +++++++ .../tdis/tracking/BuildTelescopeDetector.hpp | 57 +++ source/tdis/tracking/CKFTrackingFunction.cc | 117 ++++--- .../tdis/tracking/ReconstructedHitFactory.h | 91 ++++- source/tdis/tracking/TelescopeDetector.cpp | 62 ++++ source/tdis/tracking/TelescopeDetector.hpp | 64 ++++ .../tracking/TelescopeDetectorElement.cpp | 42 +++ .../tracking/TelescopeDetectorElement.hpp | 162 +++++++++ 16 files changed, 1123 insertions(+), 106 deletions(-) create mode 100644 source/tdis/tracking/BuildTelescopeDetector.cpp create mode 100644 source/tdis/tracking/BuildTelescopeDetector.hpp create mode 100644 source/tdis/tracking/TelescopeDetector.cpp create mode 100644 source/tdis/tracking/TelescopeDetector.hpp create mode 100644 source/tdis/tracking/TelescopeDetectorElement.cpp create mode 100644 source/tdis/tracking/TelescopeDetectorElement.hpp diff --git a/docker/tdis-pre/Dockerfile b/docker/tdis-pre/Dockerfile index 0a698f4..84cd79e 100644 --- a/docker/tdis-pre/Dockerfile +++ b/docker/tdis-pre/Dockerfile @@ -14,7 +14,7 @@ ARG BUILD_THREADS=8 # Software versions ARG VERSION_CERN_ROOT=v6-32-06 -ARG VERSION_ACTS=v37.2.0 +ARG VERSION_ACTS=v37.4.0 ARG VERSION_PODIO=v01-01 ARG VERSION_JANA2=v2.3.2 @@ -93,7 +93,7 @@ RUN edpm config root branch=${VERSION_CERN_ROOT} &&\ RUN edpm config acts branch=${VERSION_ACTS} &&\ - edpm config acts cmake_flags='-DACTS_BUILD_PLUGIN_TGEO=ON -DACTS_BUILD_PLUGIN_DD4HEP=OFF -DACTS_BUILD_PLUGIN_JSON=ON -DACTS_BUILD_PLUGIN_ACTSVG=OFF' &&\ + edpm config acts cmake_flags='-DACTS_BUILD_EXAMPLES_PYTHON_BINDINGS=ON -DACTS_BUILD_PLUGIN_TGEO=ON -DACTS_BUILD_PLUGIN_DD4HEP=OFF -DACTS_BUILD_PLUGIN_JSON=ON -DACTS_BUILD_PLUGIN_ACTSVG=OFF' &&\ edpm install acts RUN edpm config podio branch=${VERSION_PODIO} &&\ diff --git a/source/tdis/CKFTracking.h b/source/tdis/CKFTracking.h index 6373488..df67aec 100644 --- a/source/tdis/CKFTracking.h +++ b/source/tdis/CKFTracking.h @@ -38,13 +38,7 @@ namespace eicrecon { class CKFTracking: public WithPodConfig { public: - /// Track finder function that takes input measurements, initial trackstate - /// and track finder options and returns some track-finder-specific result. - using TrackFinderOptions = - Acts::CombinatorialKalmanFilterOptions; - using TrackFinderResult = - Acts::Result>; + /// Find function that takes the above parameters /// @note This is separated into a virtual interface to keep compilation units diff --git a/source/tdis/CMakeLists.txt b/source/tdis/CMakeLists.txt index 62765f0..5faa935 100644 --- a/source/tdis/CMakeLists.txt +++ b/source/tdis/CMakeLists.txt @@ -45,6 +45,14 @@ add_executable(tdis tracking/ActsGeometryService.cc tracking/ActsGeometryService.h tracking/ReconstructedHitFactory.h + + tracking/BuildTelescopeDetector.cpp + tracking/BuildTelescopeDetector.hpp + tracking/TelescopeDetector.cpp + tracking/TelescopeDetector.hpp + tracking/TelescopeDetectorElement.cpp + tracking/TelescopeDetectorElement.hpp + # tracking/CKFTrackingFunction.cc # tracking/DD4hepBField.h # /tracking/DD4hepBField.cc ) @@ -73,7 +81,39 @@ set_target_properties(tdis PROPERTIES INSTALL_RPATH_USE_LINK_PATH TRUE) target_include_directories(tdis PRIVATE "podio_model") target_include_directories(tdis SYSTEM PUBLIC ${JANA_INCLUDE_DIR} ${ROOT_INCLUDE_DIRS}) +# ----------- Configure ACTS ExamplesLibrary -------- +# ExamplesLibrary actually creates ACTS event model +# Get ActsExamples base +get_target_property(ActsCore_LOCATION ActsCore LOCATION) +get_filename_component(ActsCore_PATH ${ActsCore_LOCATION} DIRECTORY) +set(ActsExamples_LIB ${ActsCore_PATH}/${CMAKE_SHARED_LIBRARY_PREFIX}ActsExamplesFramework${CMAKE_SHARED_LIBRARY_SUFFIX}) + +# If acts is installed to / we are good but if acts is installed in some directory with structure /bin /include /lib +# we need to get this directory and try to include include +get_filename_component(Acts_HOME "${ActsCore_PATH}" DIRECTORY) +set(Acts_HOME_INCLUDE "${Acts_HOME}/include") + +# List all ACTS variables +message(STATUS "ACTS List all variables : ${ActsExample_LIB}") +get_cmake_property(_variableNames VARIABLES) +foreach (_variableName ${_variableNames}) + if (_variableName MATCHES "^Acts") + message(STATUS " ${_variableName} = ${${_variableName}}") + endif() +endforeach() + +# Add examples library and includes +# Check if the directory exists and is accessible on the disk +if(EXISTS "${Acts_HOME_INCLUDE}") + # Add the directory to the target's include directories only if it exists + target_include_directories(tdis PUBLIC SYSTEM "${Acts_HOME_INCLUDE}") + message(STATUS "Added include directory: ${Acts_HOME_INCLUDE}") +else() + message(STATUS "Directory does not exist or is not accessible: ${Acts_HOME_INCLUDE}. It is OK if ACTS is installed in /") +endif() + +target_include_directories(tdis PUBLIC ${CMAKE_CURRENT_LIST_DIR} "${CMAKE_CURRENT_LIST_DIR}/..") # ----------- install destination ------------------- install(TARGETS tdis DESTINATION bin) diff --git a/source/tdis/io/DigitizedDataEventSource.hpp b/source/tdis/io/DigitizedDataEventSource.hpp index 589ab40..d810a7e 100644 --- a/source/tdis/io/DigitizedDataEventSource.hpp +++ b/source/tdis/io/DigitizedDataEventSource.hpp @@ -59,12 +59,15 @@ namespace tdis::io { /** POD structure for readout hits **/ struct DigitizedReadoutHit { - double time; // - Time of arrival at Pad (ns) - double adc; // - Amplitude (ADC bin of sample) - int ring; // - Ring (id of rin, 0 is innermost). - int pad; // - Pad (id of pad, 0 is at or closest to phi=0 and numbering is clockwise). - int plane; // - Plane(id of z plane from 0 upstream to 9 downstream) - double zToGem; // - ZtoGEM (m) + double time; // - Time of arrival at Pad (ns) + double adc; // - Amplitude (ADC bin of sample) + int ring; // - Ring (id of rin, 0 is innermost). + int pad; // - Pad (id of pad, 0 is at or closest to phi=0 and numbering is clockwise). + int plane; // - Plane(id of z plane from 0 upstream to 9 downstream) + double zToGem; // - ZtoGEM (m) + double true_x; // - True hit x info (quiet_NaN() if not provided) + double true_y; // - True hit y info (quiet_NaN() if not provided) + double true_z; // - True hit z info (quiet_NaN() if not provided) }; /** POD structure for readout track **/ @@ -240,12 +243,31 @@ namespace tdis::io { return false; } - result.time = std::stod(tokens[0]); // - Time of arrival at Pad (ns) - result.adc = std::stod(tokens[1]); // - Amplitude (ADC bin of sample) - result.ring = std::stoi(tokens[2]); // - Ring (id of rin, 0 is innermost). - result.pad = std::stoi(tokens[3]); // - Pad (id of pad, 0 is at or closest to phi=0 and numbering is clockwise). - result.plane = std::stoi(tokens[4]); // - Plane(id of z plane from 0 upstream to 9 downstream) - result.zToGem = std::stod(tokens[5]); // - ZtoGEM (m) + if(tokens.size() == 6) { + // Files with no true X Y Z hit info + result.time = std::stod(tokens[0]); // - Time of arrival at Pad (ns) + result.adc = std::stod(tokens[1]); // - Amplitude (ADC bin of sample) + result.ring = std::stoi(tokens[2]); // - Ring (id of rin, 0 is innermost). + result.pad = std::stoi(tokens[3]); // - Pad (id of pad, 0 is at or closest to phi=0 and numbering is clockwise). + result.plane = std::stoi(tokens[4]); // - Plane(id of z plane from 0 upstream to 9 downstream) + result.zToGem = std::stod(tokens[5]); // - ZtoGEM (m) + + // True hit information is not set + result.true_x = std::numeric_limits::quiet_NaN(); + result.true_y = std::numeric_limits::quiet_NaN(); + result.true_z = std::numeric_limits::quiet_NaN(); + } else { + result.time = std::stod(tokens[0]); // - Time of arrival at Pad (ns) + result.adc = std::stod(tokens[1]); // - Amplitude (ADC bin of sample) + result.true_x = std::stod(tokens[2]); // True X Y Z of hit + result.true_y = std::stod(tokens[3]); + result.true_z = std::stod(tokens[4]); + result.ring = std::stoi(tokens[5]); // - Ring (id of rin, 0 is innermost). + result.pad = std::stoi(tokens[6]); // - Pad (id of pad, 0 is at or closest to phi=0 and numbering is clockwise). + result.plane = std::stoi(tokens[7]); // - Plane(id of z plane from 0 upstream to 9 downstream) + result.zToGem = std::stod(tokens[8]); // - ZtoGEM (m) + } + return true; } diff --git a/source/tdis/layout.yaml b/source/tdis/layout.yaml index a63ca6f..4c2c491 100644 --- a/source/tdis/layout.yaml +++ b/source/tdis/layout.yaml @@ -289,6 +289,9 @@ datatypes : - int pad // - Pad (id of pad, 0 is at or closest to phi=0 and numbering is clockwise). - int plane // - Plane(id of z plane from 0 upstream to 9 downstream) - double zToGem // - ZtoGEM (m) + - double true_x // - True hit x info (quiet_NaN() if not provided) + - double true_y // - True hit y info (quiet_NaN() if not provided) + - double true_z // - True hit z info (quiet_NaN() if not provided) edm4eic::TrackerHit: Description: "Tracker hit (reconstructed from Raw)" diff --git a/source/tdis/tdis_main.cpp b/source/tdis/tdis_main.cpp index 6fc9c81..8954c66 100644 --- a/source/tdis/tdis_main.cpp +++ b/source/tdis/tdis_main.cpp @@ -15,9 +15,11 @@ #include #include "CLI/CLI.hpp" -#include "io/PodioWriteProcessor.hpp" #include "io/DigitizedDataEventSource.hpp" +#include "io/PodioWriteProcessor.hpp" #include "services/LogService.hpp" +#include "tracking/ActsGeometryService.h" +#include "tracking/ReconstructedHitFactory.h" struct ProgramArguments { std::map params; @@ -106,12 +108,15 @@ int main(int argc, char* argv[]) { parameterManager->SetParameter(name, value); } + JApplication app(parameterManager); + // Register services: + app.ProvideService(std::make_shared(&app)); + app.ProvideService(std::make_shared()); - JApplication app(parameterManager); app.Add(new JEventSourceGeneratorT); app.Add(new tdis::io::PodioWriteProcessor(&app)); - app.ProvideService(std::make_shared(&app)); + app.Add(new JFactoryGeneratorT); // app.Add(new JEventProcessorPodio); // app.Add(new JFactoryGeneratorT()); // app.Add(new JFactoryGeneratorT()); diff --git a/source/tdis/tracking/ActsGeometryService.cc b/source/tdis/tracking/ActsGeometryService.cc index 3ae1ad2..4f21c57 100644 --- a/source/tdis/tracking/ActsGeometryService.cc +++ b/source/tdis/tracking/ActsGeometryService.cc @@ -7,17 +7,24 @@ #include #include +#include #include #include +#include #include #include +#include +#include #include #include #include #include +#include #include "ActsGeometryProvider.h" +#include "BuildTelescopeDetector.hpp" +#include "TelescopeDetector.hpp" #include "services/LogService.hpp" // Formatter for Eigen matrices @@ -38,7 +45,217 @@ struct fmt::formatter< #endif // FMT_VERSION >= 90000 -void ActsGeometryService::Init() { +namespace { + // Function to recursively find the node by volume name + inline TGeoNode* findNodeRecursive(TGeoNode* currentNode, const char* volumeName) { + + // Check if the current node's volume matches the name + if (std::strcmp(currentNode->GetVolume()->GetName(), volumeName) == 0) { + return currentNode; + } + // Recursively search in the daughter nodes + int nDaughters = currentNode->GetNdaughters(); + for (int i = 0; i < nDaughters; ++i) { + TGeoNode* daughterNode = currentNode->GetDaughter(i); + TGeoNode* foundNode = findNodeRecursive(daughterNode, volumeName); + if (foundNode != nullptr) { + return foundNode; + } + } + return nullptr; // Not found in this branch + } + + + + void printNodeTree(TGeoNode* currentNode, bool printVolumes = false, int level = 0) { + // Print spaces corresponding to the level + for (int i = 0; i < level; ++i) { + std::cout << " "; // Two spaces per level + } + + // Print the node's name or volume's name + if (printVolumes) { + std::cout << currentNode->GetVolume()->GetName() << std::endl; + } else { + std::cout << currentNode->GetName() << std::endl; + } + + // Recursively print daughter nodes + int nDaughters = currentNode->GetNdaughters(); + for (int i = 0; i < nDaughters; ++i) { + TGeoNode* daughterNode = currentNode->GetDaughter(i); + printNodeTree(daughterNode, printVolumes, level + 1); + } + } + + /** + * + */ + void findNodesWithPrefix(TGeoNode* currentNode, const std::string& prefix, std::vector& results, bool searchVolumes = false) { + std::string_view name = searchVolumes ? currentNode->GetVolume()->GetName() + : currentNode->GetName(); + + if (name.starts_with(prefix)) { + results.push_back(currentNode); + } + + // Recursively search in the daughter nodes + int nDaughters = currentNode->GetNdaughters(); + for (int i = 0; i < nDaughters; ++i) { + TGeoNode* daughterNode = currentNode->GetDaughter(i); + findNodesWithPrefix(daughterNode, prefix, results, searchVolumes); + } + } + + double roundTo2DecimalPlaces(double value) { + return std::round(value * 100.0) / 100.0; + } + + +#include +#include +#include "TGeoManager.h" +#include "TGeoNode.h" +#include "TGeoVolume.h" +#include "TGeoShape.h" +#include "TGeoMatrix.h" + + // Recursive function to find the node and compute its global matrix + bool findNodeGlobalMatrix(TGeoNode* currentNode, TGeoNode* targetNode, TGeoHMatrix& globalMatrix) { + if (!currentNode) { + return false; + } + + // Save the current global matrix + TGeoHMatrix savedMatrix = globalMatrix; + + // Get the local transformation matrix of the current node + TGeoMatrix* localMatrix = currentNode->GetMatrix(); + if (localMatrix) { + // Multiply the global matrix by the local matrix + globalMatrix.Multiply(localMatrix); + } + + // Check if currentNode is the targetNode + if (currentNode == targetNode) { + // Found the node + return true; + } + + // Recursively search in daughter nodes + int nDaughters = currentNode->GetNdaughters(); + for (int i = 0; i < nDaughters; ++i) { + TGeoNode* daughterNode = currentNode->GetDaughter(i); + + // Create a copy of the current global matrix for the daughter + TGeoHMatrix daughterGlobalMatrix = globalMatrix; + + if (findNodeGlobalMatrix(daughterNode, targetNode, daughterGlobalMatrix)) { + // If the node is found in this branch, update the global matrix + globalMatrix = daughterGlobalMatrix; + return true; + } + } + + // Node not found in this branch, restore the global matrix + globalMatrix = savedMatrix; + return false; + } + + std::tuple getGlobalPosition(TGeoNode* node) { + // Initialize the global matrix as identity + TGeoHMatrix globalMatrix; + + // Start from the top node + TGeoNode* topNode = gGeoManager->GetTopNode(); + + + // Find the node and compute its global matrix + if (findNodeGlobalMatrix(topNode, node, globalMatrix)) { + // Extract the translation (position) from the global matrix + const Double_t* translation = globalMatrix.GetTranslation(); + + // Print the global position + return {translation[0], translation[1], translation[2]}; + } + + auto message = fmt::format("Node {} not found in the geometry tree.", node->GetName()); + throw std::runtime_error(message); + } + + void printNodeGlobalPosition(TGeoNode* node) { + if (!node) { + std::cerr << "Invalid node provided." << std::endl; + return; + } + + // Initialize the global matrix as identity + TGeoHMatrix globalMatrix; + + // Start from the top node + TGeoNode* topNode = gGeoManager->GetTopNode(); + if (!topNode) { + std::cerr << "Top node not found in the geometry manager." << std::endl; + return; + } + + // Find the node and compute its global matrix + if (findNodeGlobalMatrix(topNode, node, globalMatrix)) { + // Extract the translation (position) from the global matrix + const Double_t* translation = globalMatrix.GetTranslation(); + + // Print the global position + std::cout << "Global position of the node '" << node->GetName() << "': (" + << translation[0] << ", " + << translation[1] << ", " + << translation[2] << ")" << std::endl; + } else { + std::cerr << "Node not found in the geometry tree." << std::endl; + } + } + + void printNodeInfo(TGeoNode* node) { + if (!node) { + std::cerr << "Invalid node provided." << std::endl; + return; + } + + // Get the navigator + TGeoNavigator* navigator = gGeoManager->GetCurrentNavigator(); + if (!navigator) { + // Create a navigator if one doesn't exist + navigator = gGeoManager->AddNavigator(); + } + + // Save the current position + TString currentPath = navigator->GetPath(); + + // Get the volume and shape associated with the node + TGeoVolume* volume = node->GetVolume(); + TGeoShape* shape = volume->GetShape(); + + // Get and print the shape type + std::string shapeType = shape->ClassName(); + std::cout << "Shape Type: " << shapeType << std::endl; + + // Retrieve and print dimensions based on the shape type + if (shapeType == "TGeoTube" || shapeType == "TGeoTubeSeg") { + TGeoTube* tube = dynamic_cast(shape); + if (tube) { + std::cout << "Inner Radius (Rmin): " << tube->GetRmin() << std::endl; + std::cout << "Outer Radius (Rmax): " << tube->GetRmax() << std::endl; + std::cout << "Half Length Z (Dz): " << tube->GetDz() << std::endl; + } + } + // Handle other shape types as needed + + // Restore the previous position + navigator->cd(currentPath.Data()); + } +} + + +void tdis::tracking::ActsGeometryService::Init() { m_log = m_service_log->logger("acts"); m_init_log = m_service_log->logger("acts_init"); @@ -68,13 +285,119 @@ void ActsGeometryService::Init() { // Set up the converter first Acts::MaterialMapJsonConverter::Config jsonGeoConvConfig; // Set up the json-based decorator - materialDeco = std::make_shared(jsonGeoConvConfig, m_material_map_file(), acts_init_log_level); + try { + materialDeco = std::make_shared(jsonGeoConvConfig, m_material_map_file(), acts_init_log_level); + } + catch (const std::exception& e) { + m_init_log->error("Failed to load materials map: {}", e.what()); + exit(1); // TODO this is due to JANA2 issue #381. Remove after is fixed + } } // Convert DD4hep geometry to ACTS m_init_log->info("Converting TGeo geometry to ACTS..."); auto logger = eicrecon::getSpdlogLogger("CONV", m_log); + m_log->info("Loading geometry file: "); + m_log->info(" '{}'", m_tgeo_file()); + + m_tgeo_manager = TGeoManager::Import(m_tgeo_file().c_str()); + if(!m_tgeo_manager) { + m_init_log->error("Failed to load GeometryFile: TGeoManager::Import returned NULL"); + exit(1); // TODO this is due to JANA2 issue #381. Remove after is fixed + } + + if(!m_tgeo_manager->GetTopNode()) { + m_init_log->error("Failed to load GeometryFile: TGeoManager::GetTopNode returned NULL"); + exit(1); // TODO this is due to JANA2 issue #381. Remove after is fixed + } + + if(!m_tgeo_manager->GetTopVolume()) { + m_init_log->error("Failed to load GeometryFile: TGeoManager::GetTopVolume returned NULL"); + exit(1); // TODO this is due to JANA2 issue #381. Remove after is fixed + } + + printNodeTree(m_tgeo_manager->GetTopNode(), /*printVolumes*/ false); + + std::vector disk_nodes; + + findNodesWithPrefix(m_tgeo_manager->GetTopNode(), "mTPCReadoutDisc", disk_nodes); + + + + std::vector positions; + std::vector stereos; + std::array offsets{{0, 0}}; + std::array bounds{{25, 100}}; + double thickness{0.1}; + int surfaceType{(int) ActsExamples::Telescope::TelescopeSurfaceType::Disc}; + + + m_init_log->info("Found readout disks: {}", disk_nodes.size()); + for (auto node : disk_nodes) { + + auto [x,y,z] = getGlobalPosition(node); + // Get the volume and shape associated with the node + TGeoVolume* volume = node->GetVolume(); + TGeoShape* shape = volume->GetShape(); + std::string shapeType = shape->ClassName(); + + // Check it is TGeoTube as expected + if (shapeType != "TGeoTube") { + m_init_log->warn("For TGeoNode('{}') TGeoVolume ('{}') the shape not a TGeoTube as expected but: ", node->GetName(), volume->GetName(), shapeType); + m_init_log->warn("We consider this as recoverable, but it might be if all readout disks are defined with other shape or something. Check this!"); + continue; + } + + TGeoTube* tube = dynamic_cast(shape); + double r_min = tube->GetRmin(); + double r_max = tube->GetRmax(); + double dz = tube->GetDz(); + + // We just copy it, considering the disks are the same + bounds[0] = r_min * Acts::UnitConstants::cm; + bounds[1] = r_max * Acts::UnitConstants::cm; + thickness = dz * Acts::UnitConstants::cm * 0.1; + + m_log->info(fmt::format(" Position: ({:.4f}, {:.4f}, {:>8.4f}) cm, rmin: {:.3f} cm, rmax: {:.3f} cm, dz: {:.4f} cm, Node: '{}', Volume: '{}'", + x, y, z, r_min, r_max, dz, node->GetName(), volume->GetName())); + + positions.push_back(z * Acts::UnitConstants::cm); + stereos.push_back(0); + } + + Acts::ObjVisualization3D objVis; + + m_init_log->info("Building ACTS Geometry"); + + + + /// Return the telescope detector + std::shared_ptr gGeometry = + ActsExamples::Telescope::buildDetector( + nominalContext, // gctx is the detector element dependent geometry context + detectorStore, // detectorStore is the store for the detector element + positions, // positions are the positions of different layers in the longitudinal direction + stereos, // stereoAngles are the stereo angles of different layers, which are rotation angles around the longitudinal (normal) direction + offsets, // is the offset (u, v) of the layers in the transverse plane + bounds, // bounds is the surface bound values, i.e. halfX and halfY if plane surface, and minR and maxR if disc surface + thickness, // thickness is the material thickness of each layer + ActsExamples::Telescope::TelescopeSurfaceType::Disc, // surfaceType is the detector surface type + Acts::BinningValue::binZ); + + m_log->info("ACTS Detector Elements"); + for(auto &element: detectorStore) { + Acts::GeometryView3D::drawSurface(objVis, element->surface(), m_geometry_context, element->transform(m_geometry_context)); + } + + if(!m_obj_output_file().empty()) { + m_log->info("ACTS Exporting geometry as view to: {}", m_obj_output_file()); + objVis.write(m_obj_output_file()); + } else { + m_log->info("ACTS geometry filename defined by: '{}' is empty. NOT exporting geometry to .obj or whatever ", m_obj_output_file.m_name); + } + + exit(0); // Set ticker back diff --git a/source/tdis/tracking/ActsGeometryService.h b/source/tdis/tracking/ActsGeometryService.h index 2ca92f4..d646766 100644 --- a/source/tdis/tracking/ActsGeometryService.h +++ b/source/tdis/tracking/ActsGeometryService.h @@ -6,42 +6,66 @@ #include #include #include +#include #include +#include "TelescopeDetectorElement.hpp" + #include #include #include "ActsGeometryProvider.h" #include "services/LogService.hpp" -class ActsGeometryService : public JService -{ -public: - explicit ActsGeometryService() : JService() {} - ~ActsGeometryService() override = default; +namespace tdis::tracking { + class ActsGeometryService : public JService + { + public: + explicit ActsGeometryService() : JService() {} + ~ActsGeometryService() override = default; + + + void Init() override; + + TGeoManager* GetGeoManager() const { return m_tgeo_manager; } + + TGeoVolume * GetTopVolume() const { return m_tgeo_manager->GetTopVolume(); } + + TGeoNode * GetTopNode() const { return m_tgeo_manager->GetTopNode(); } + + + Acts::GeometryContext& GetActsGeometryContext() { + return m_geometry_context; + } + + + private: - void Init() override; + Parameter m_tgeo_file {this, "acts:geometry", "g4sbs_mtpc.root", "TGeo filename with geometry for ACTS"}; + Parameter m_material_map_file {this, "acts:material_map", "", "JSON/CBOR material map file path"}; + Parameter m_obj_output_file {this, "acts:output_obj", "", "Output file name to dump ACTS converted geometry as OBJ"}; + Parameter m_ply_output_file {this, "acts:output_ply", "", "Output file name to dump ACTS converted geometry as PLY"}; -protected: + Service m_service_log {this}; -private: + // General acts log + std::shared_ptr m_log; + /// Logger that is used for geometry initialization + /// By default its level the same as ACTS general logger (m_log) + /// But it might be customized to solely printout geometry information + std::shared_ptr m_init_log; - Parameter m_output_filename {this, "acts:geometry", "g4sbs_mtpc.root", "TGeo filename with geometry for ACTS"}; - Parameter m_material_map_file {this, "acts:material_map", "material_map.json", "JSON/CBOR material map file path"}; - Parameter m_obj_output_file {this, "acts:output_obj", "", "Output file name to dump ACTS converted geometry as OBJ"}; - Parameter m_ply_output_file {this, "acts:output_ply", "", "Output file name to dump ACTS converted geometry as PLY"}; + /// Root TGeo Manater for TGeo Geometry + TGeoManager* m_tgeo_manager; - Service m_service_log {this}; + Acts::GeometryContext m_geometry_context = Acts::GeometryContext(); + ActsExamples::Telescope::TelescopeDetectorElement::ContextType nominalContext; - // General acts log - std::shared_ptr m_log; + std::vector> detectorStore; - /// Logger that is used for geometry initialization - /// By default its level the same as ACTS general logger (m_log) - /// But it might be customized to solely printout geometry information - std::shared_ptr m_init_log; -}; + }; +} // namespace tdis::tracking \ No newline at end of file diff --git a/source/tdis/tracking/BuildTelescopeDetector.cpp b/source/tdis/tracking/BuildTelescopeDetector.cpp new file mode 100644 index 0000000..c237127 --- /dev/null +++ b/source/tdis/tracking/BuildTelescopeDetector.cpp @@ -0,0 +1,131 @@ +// This file is part of the Acts project. +// +// Copyright (C) 2020-2021 CERN for the benefit of the Acts project +// +// This Source Code Form is subject to the terms of the Mozilla Public +// License, v. 2.0. If a copy of the MPL was not distributed with this +// file, You can obtain one at http://mozilla.org/MPL/2.0/. + +#include "BuildTelescopeDetector.hpp" + +#include "Acts/Definitions/Algebra.hpp" +#include "Acts/Definitions/Units.hpp" +#include "Acts/Geometry/CuboidVolumeBounds.hpp" +#include "Acts/Geometry/CylinderVolumeBounds.hpp" +#include "Acts/Geometry/DiscLayer.hpp" +#include "Acts/Geometry/GeometryContext.hpp" +#include "Acts/Geometry/ILayerArrayCreator.hpp" +#include "Acts/Geometry/ITrackingVolumeHelper.hpp" +#include "Acts/Geometry/LayerArrayCreator.hpp" +#include "Acts/Geometry/PlaneLayer.hpp" +#include "Acts/Geometry/TrackingGeometry.hpp" +#include "Acts/Geometry/TrackingVolume.hpp" +#include "Acts/Material/HomogeneousSurfaceMaterial.hpp" +#include "Acts/Material/Material.hpp" +#include "Acts/Material/MaterialSlab.hpp" +#include "Acts/Surfaces/RadialBounds.hpp" +#include "Acts/Surfaces/RectangleBounds.hpp" +#include "Acts/Surfaces/Surface.hpp" +#include "Acts/Surfaces/SurfaceArray.hpp" +#include "Acts/Utilities/Logger.hpp" + +#include +#include +#include + +std::unique_ptr +ActsExamples::Telescope::buildDetector(const typename ActsExamples::Telescope::TelescopeDetectorElement::ContextType &gctx, std::vector > &detectorStore, const std::vector &positions, const std::vector &stereoAngles, const std::array &offsets, const std::array &bounds, double thickness, ActsExamples::Telescope::TelescopeSurfaceType surfaceType, Acts::BinningValue binValue) { + using namespace Acts::UnitLiterals; + + // The rectangle bounds for plane surface + const auto pBounds = std::make_shared(bounds[0], bounds[1]); + // The radial bounds for disc surface + const auto rBounds = std::make_shared(bounds[0], bounds[1]); + + // Material of the surfaces + Acts::Material silicon = Acts::Material::fromMassDensity(9.370_cm, 46.52_cm, 28.0855, 14, 2.329_g / 1_cm3); + Acts::MaterialSlab matProp(silicon, thickness); + const auto surfaceMaterial = std::make_shared(matProp); + + // Construct the rotation + // This assumes the binValue is binX, binY or binZ. No reset is necessary in + // case of binZ + Acts::RotationMatrix3 rotation = Acts::RotationMatrix3::Identity(); + if (binValue == Acts::BinningValue::binX) { + rotation.col(0) = Acts::Vector3(0, 0, -1); + rotation.col(1) = Acts::Vector3(0, 1, 0); + rotation.col(2) = Acts::Vector3(1, 0, 0); + } else if (binValue == Acts::BinningValue::binY) { + rotation.col(0) = Acts::Vector3(1, 0, 0); + rotation.col(1) = Acts::Vector3(0, 0, -1); + rotation.col(2) = Acts::Vector3(0, 1, 0); + } + + // Construct the surfaces and layers + std::size_t nLayers = positions.size(); + std::vector layers(nLayers); + for (unsigned int i = 0; i < nLayers; ++i) { + // The translation without rotation yet + Acts::Translation3 trans(offsets[0], offsets[1], positions[i]); + // The entire transformation (the coordinate system, whose center is defined + // by trans, will be rotated as well) + Acts::Transform3 trafo(rotation * trans); + + // rotate around local z axis by stereo angle + auto stereo = stereoAngles[i]; + trafo *= Acts::AngleAxis3(stereo, Acts::Vector3::UnitZ()); + + // Create the detector element + std::shared_ptr detElement = nullptr; + if (surfaceType == TelescopeSurfaceType::Plane) { + detElement = std::make_shared(std::make_shared(trafo), pBounds, 1._um, surfaceMaterial); + } else { + detElement = std::make_shared(std::make_shared(trafo), rBounds, 1._um, surfaceMaterial); + } + // Get the surface + auto surface = detElement->surface().getSharedPtr(); + // Add the detector element to the detector store + detectorStore.push_back(std::move(detElement)); + // Construct the surface array (one surface contained) + std::unique_ptr surArray(new Acts::SurfaceArray(surface)); + // Construct the layer + if (surfaceType == TelescopeSurfaceType::Plane) { + layers[i] = Acts::PlaneLayer::create(trafo, pBounds, std::move(surArray), 1._mm); + } else { + layers[i] = Acts::DiscLayer::create(trafo, rBounds, std::move(surArray), 1._mm); + } + // Associate the layer to the surface + auto mutableSurface = const_cast(surface.get()); + mutableSurface->associateLayer(*layers[i]); + } + + // The volume transform + Acts::Translation3 transVol(offsets[0], offsets[1], (positions.front() + positions.back()) * 0.5); + Acts::Transform3 trafoVol(rotation * transVol); + + // The volume bounds is set to be a bit larger than either cubic with planes + // or cylinder with discs + auto length = positions.back() - positions.front(); + std::shared_ptr boundsVol = nullptr; + if (surfaceType == TelescopeSurfaceType::Plane) { + boundsVol = std::make_shared(bounds[0] + 5._mm, bounds[1] + 5._mm, length + 10._mm); + } else { + boundsVol = std::make_shared(std::max(bounds[0] - 5.0_mm, 0.), bounds[1] + 5._mm, length + 10._mm); + } + + Acts::LayerArrayCreator::Config lacConfig; + Acts::LayerArrayCreator layArrCreator(lacConfig, Acts::getDefaultLogger("LayerArrayCreator", Acts::Logging::INFO)); + Acts::LayerVector layVec; + for (unsigned int i = 0; i < nLayers; i++) { + layVec.push_back(layers[i]); + } + // Create the layer array + Acts::GeometryContext genGctx{gctx}; + std::unique_ptr layArr(layArrCreator.layerArray(genGctx, layVec, positions.front() - 2._mm, positions.back() + 2._mm, Acts::BinningType::arbitrary, binValue)); + + // Build the tracking volume + auto trackVolume = std::make_shared(trafoVol, boundsVol, nullptr, std::move(layArr), nullptr, Acts::MutableTrackingVolumeVector{}, "Telescope"); + + // Build and return tracking geometry + return std::make_unique(trackVolume); +} diff --git a/source/tdis/tracking/BuildTelescopeDetector.hpp b/source/tdis/tracking/BuildTelescopeDetector.hpp new file mode 100644 index 0000000..f64c5c3 --- /dev/null +++ b/source/tdis/tracking/BuildTelescopeDetector.hpp @@ -0,0 +1,57 @@ +// This file is part of the Acts project. +// +// Copyright (C) 2020-2021 CERN for the benefit of the Acts project +// +// This Source Code Form is subject to the terms of the Mozilla Public +// License, v. 2.0. If a copy of the MPL was not distributed with this +// file, You can obtain one at http://mozilla.org/MPL/2.0/. + +#pragma once + +#include "Acts/Geometry/TrackingGeometry.hpp" +#include "Acts/Utilities/BinUtility.hpp" +#include "Acts/Utilities/BinningType.hpp" +#include "TelescopeDetectorElement.hpp" + +#include +#include +#include + +namespace Acts { +class TrackingGeometry; +} // namespace Acts + +namespace ActsExamples::Telescope { + +/// The telescope detector surface type +enum class TelescopeSurfaceType { + Plane = 0, + Disc = 1, +}; + +/// Global method to build the telescope tracking geometry +/// +/// @param gctx is the detector element dependent geometry context +/// @param detectorStore is the store for the detector element +/// @param positions are the positions of different layers in the longitudinal +/// direction +/// @param stereoAngles are the stereo angles of different layers, which are +/// rotation angles around the longitudinal (normal) +/// direction +/// @param offsets is the offset (u, v) of the layers in the transverse plane +/// @param bounds is the surface bound values, i.e. halfX and halfY if plane +/// surface, and minR and maxR if disc surface +/// @param thickness is the material thickness of each layer +/// @param surfaceType is the detector surface type +/// @param binValue indicates which axis the detector surface normals are +/// parallel to +std::unique_ptr buildDetector( + const typename TelescopeDetectorElement::ContextType& gctx, + std::vector>& detectorStore, + const std::vector& positions, + const std::vector& stereoAngles, + const std::array& offsets, const std::array& bounds, + double thickness, TelescopeSurfaceType surfaceType, + Acts::BinningValue binValue = Acts::BinningValue::binZ); + +} // namespace ActsExamples::Telescope diff --git a/source/tdis/tracking/CKFTrackingFunction.cc b/source/tdis/tracking/CKFTrackingFunction.cc index 5e748ff..4a18fec 100644 --- a/source/tdis/tracking/CKFTrackingFunction.cc +++ b/source/tdis/tracking/CKFTrackingFunction.cc @@ -20,71 +20,88 @@ #include #include #include -#include #include #include +#include #include #include #include "ActsExamples/EventData/Track.hpp" -#include "CKFTracking.h" +// #include "CKFTracking.h" -namespace eicrecon{ +namespace tdis { - using Updater = Acts::GainMatrixUpdater; - using Smoother = Acts::GainMatrixSmoother; - using Stepper = Acts::EigenStepper<>; - using Navigator = Acts::Navigator; - using Propagator = Acts::Propagator; - using CKF = - Acts::CombinatorialKalmanFilter; + using Updater = Acts::GainMatrixUpdater; + using Smoother = Acts::GainMatrixSmoother; - using TrackContainer = - Acts::TrackContainer; + using Stepper = Acts::EigenStepper<>; + using Navigator = Acts::Navigator; + using Propagator = Acts::Propagator; - /** Finder implementation . - * - * \ingroup track - */ - struct CKFTrackingFunctionImpl - : public eicrecon::CKFTracking::CKFTrackingFunction { - CKF trackFinder; + using CombinatorialKalmanFilter = Acts::CombinatorialKalmanFilter; - CKFTrackingFunctionImpl(CKF&& f) : trackFinder(std::move(f)) {} - - eicrecon::CKFTracking::TrackFinderResult operator()( - const ActsExamples::TrackParameters& initialParameters, - const eicrecon::CKFTracking::TrackFinderOptions& options, - TrackContainer& tracks) const override { - return trackFinder.findTracks(initialParameters, options, tracks); - }; - }; + using TrackContainer = Acts::TrackContainer; -} // namespace + /// Track finder function that takes input measurements, initial track state + /// and track finder options and returns some track-finder-specific result. + using TrackFinderOptions = Acts::CombinatorialKalmanFilterOptions; + using TrackFinderResult = Acts::Result>; -namespace eicrecon { + /// Find function that takes the above parameters + /// @note This is separated into a virtual interface to keep compilation units + /// small + class CKFTrackingFunction { + public: + virtual ~CKFTrackingFunction() = default; - std::shared_ptr - CKFTracking::makeCKFTrackingFunction( - std::shared_ptr trackingGeometry, - std::shared_ptr magneticField, - const Acts::Logger& logger) - { - Stepper stepper(std::move(magneticField)); - Navigator::Config cfg{trackingGeometry}; - cfg.resolvePassive = false; - cfg.resolveMaterial = true; - cfg.resolveSensitive = true; - Navigator navigator(cfg); - - Propagator propagator(std::move(stepper), std::move(navigator)); - CKF trackFinder(std::move(propagator), logger.cloneWithSuffix("CKF")); + virtual TrackFinderResult operator()(const ActsExamples::TrackParameters&, const TrackFinderOptions&, ActsExamples::TrackContainer&) const = 0; + }; - // build the track finder functions. owns the track finder object. - return std::make_shared(std::move(trackFinder)); - } + /// Create the track finder function implementation. + /// The magnetic field is intentionally given by-value since the variantresults + /// contains shared_ptr anyways. + static std::shared_ptr makeCKFTrackingFunction( + std::shared_ptr trackingGeometry, + std::shared_ptr magneticField, + const Acts::Logger& logger); + + /** Finder implementation . + * + * \ingroup track + */ + struct CKFTrackingFunctionImpl : public tdis::CKFTrackingFunction { + CombinatorialKalmanFilter trackFinder; + + CKFTrackingFunctionImpl(CombinatorialKalmanFilter&& f) : trackFinder(std::move(f)) {} + + tdis::TrackFinderResult operator()( + const ActsExamples::TrackParameters& initialParameters, + const tdis::TrackFinderOptions& options, + TrackContainer& tracks) const override { + return trackFinder.findTracks(initialParameters, options, tracks); + }; + }; -} // namespace eicrecon::Reco + std::shared_ptr makeCKFTrackingFunction( + std::shared_ptr trackingGeometry, + std::shared_ptr magneticField, + const Acts::Logger& logger) + { + Stepper stepper(std::move(magneticField)); + Navigator::Config cfg{trackingGeometry}; + cfg.resolvePassive = false; + cfg.resolveMaterial = true; + cfg.resolveSensitive = true; + Navigator navigator(cfg); + + Propagator propagator(std::move(stepper), std::move(navigator)); + CombinatorialKalmanFilter trackFinder(std::move(propagator), logger.cloneWithSuffix("CKF")); + + // build the track finder functions. owns the track finder object. + return std::make_shared(std::move(trackFinder)); + } + +} // namespace tdis diff --git a/source/tdis/tracking/ReconstructedHitFactory.h b/source/tdis/tracking/ReconstructedHitFactory.h index f6c99a6..0fcc5e4 100644 --- a/source/tdis/tracking/ReconstructedHitFactory.h +++ b/source/tdis/tracking/ReconstructedHitFactory.h @@ -1,19 +1,28 @@ #pragma once -#include #include +#include + +#include "PadGeometryHelper.hpp" +#include "podio_model/DigitizedMtpcMcHit.h" +#include "podio_model/MutableTrackerHit.h" +#include "podio_model/TrackerHit.h" +#include "podio_model/TrackerHitCollection.h" namespace tdis::tracking { - struct MyClusterFactory : public JOmniFactory { + struct ReconstructedHitFactory : public JOmniFactory { - PodioInput m_protoclusters_in {this}; - PodioOutput m_clusters_out {this}; + PodioInput m_mc_hits_in {this}; + PodioOutput m_tracker_hits_out {this}; + Service m_service_geometry{this}; void Configure() { + m_service_geometry(); + } void ChangeRun(int32_t /*run_nr*/) { @@ -21,16 +30,78 @@ namespace tdis::tracking { void Execute(int32_t /*run_nr*/, uint64_t /*evt_nr*/) { - auto cs = std::make_unique(); + auto result_hits = std::make_unique(); + + for (auto mc_hit : *m_mc_hits_in()) { + + auto [x,y] = getPadCenter(mc_hit.ring(), mc_hit.pad()); + //#double z = getRing + + //std::uint64_t cellID, edm4hep::Vector3f position, edm4eic::CovDiag3f positionError, float time, float timeError, float edep, float edepError + //cluster.addClusters(protocluster); + //cs->push_back(cluster); + } + + m_tracker_hits_out() = std::move(result_hits); - for (auto protocluster : *m_protoclusters_in()) { - auto cluster = MutableExampleCluster(protocluster.energy() + 1000); - cluster.addClusters(protocluster); - cs->push_back(cluster); + /* + *auto rec_hits { std::make_unique() }; + + for (const auto& raw_hit : raw_hits) { + + auto id = raw_hit.getCellID(); + + // Get position and dimension + auto pos = m_converter->position(id); + auto dim = m_converter->cellDimensions(id); + + // >oO trace + if(m_log->level() == spdlog::level::trace) { + m_log->trace("position x={:.2f} y={:.2f} z={:.2f} [mm]: ", pos.x()/ mm, pos.y()/ mm, pos.z()/ mm); + m_log->trace("dimension size: {}", dim.size()); + for (size_t j = 0; j < std::size(dim); ++j) { + m_log->trace(" - dimension {:<5} size: {:.2}", j, dim[j]); } + } - m_clusters_out() = std::move(cs); + // Note about variance: + // The variance is used to obtain a diagonal covariance matrix. + // Note that the covariance matrix is written in DD4hep surface coordinates, + // *NOT* global position coordinates. This implies that: + // - XY segmentation: xx -> sigma_x, yy-> sigma_y, zz -> 0, tt -> 0 + // - XZ segmentation: xx -> sigma_x, yy-> sigma_z, zz -> 0, tt -> 0 + // - 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. + #if EDM4EIC_VERSION_MAJOR >= 7 + auto rec_hit = + #endif + rec_hits->create( + raw_hit.getCellID(), // Raw DD4hep cell ID + edm4hep::Vector3f{static_cast(pos.x() / mm), static_cast(pos.y() / mm), static_cast(pos.z() / mm)}, // mm + edm4eic::CovDiag3f{get_variance(dim[0] / mm), get_variance(dim[1] / mm), // variance (see note above) + std::size(dim) > 2 ? get_variance(dim[2] / mm) : 0.}, + static_cast((double)(raw_hit.getTimeStamp()) / 1000.0), // ns + m_cfg.timeResolution, // in ns + static_cast(raw_hit.getCharge() / 1.0e6), // Collected energy (GeV) + 0.0F); // Error on the energy +#if EDM4EIC_VERSION_MAJOR >= 7 + rec_hit.setRawHit(raw_hit); +#endif + * */ } + + private: + static inline double get_resolution(const double pixel_size) { + constexpr const double sqrt_12 = 3.4641016151; + return pixel_size / sqrt_12; + } + + static inline double get_variance(const double pixel_size) { + const double res = get_resolution(pixel_size); + return res * res; + } + }; } \ No newline at end of file diff --git a/source/tdis/tracking/TelescopeDetector.cpp b/source/tdis/tracking/TelescopeDetector.cpp new file mode 100644 index 0000000..ba4f7bf --- /dev/null +++ b/source/tdis/tracking/TelescopeDetector.cpp @@ -0,0 +1,62 @@ +// This file is part of the Acts project. +// +// Copyright (C) 2020-2021 CERN for the benefit of the Acts project +// +// This Source Code Form is subject to the terms of the Mozilla Public +// License, v. 2.0. If a copy of the MPL was not distributed with this +// file, You can obtain one at http://mozilla.org/MPL/2.0/. + +#include "TelescopeDetector.hpp" + +#include "Acts/Geometry/TrackingGeometry.hpp" +#include "Acts/Utilities/BinningType.hpp" +#include "BuildTelescopeDetector.hpp" +#include "TelescopeDetectorElement.hpp" + +#include +#include + +auto ActsExamples::Telescope::TelescopeDetector::finalize( + const Config& cfg, const std::shared_ptr& + /*mdecorator*/) -> std::pair { + DetectorElement::ContextType nominalContext; + + if (cfg.surfaceType > 1) { + throw std::invalid_argument( + "The surface type could either be 0 for plane surface or 1 for disc " + "surface."); + } + if (cfg.binValue > 2) { + throw std::invalid_argument("The axis value could only be 0, 1, or 2."); + } + // Check if the bounds values are valid + if (cfg.surfaceType == 1 && cfg.bounds[0] >= cfg.bounds[1]) { + throw std::invalid_argument( + "The minR should be smaller than the maxR for disc surface bounds."); + } + + if (cfg.positions.size() != cfg.stereos.size()) { + throw std::invalid_argument( + "The number of provided positions must match the number of " + "provided stereo angles."); + } + + config = cfg; + + // Sort the provided distances + std::vector positions = cfg.positions; + std::vector stereos = cfg.stereos; + std::sort(positions.begin(), positions.end()); + + /// Return the telescope detector + TrackingGeometryPtr gGeometry = ActsExamples::Telescope::buildDetector( + nominalContext, detectorStore, positions, stereos, cfg.offsets, + cfg.bounds, cfg.thickness, + static_cast( + cfg.surfaceType), + static_cast(cfg.binValue)); + ContextDecorators gContextDecorators = {}; + // return the pair of geometry and empty decorators + return std::make_pair( + std::move(gGeometry), std::move(gContextDecorators)); +} diff --git a/source/tdis/tracking/TelescopeDetector.hpp b/source/tdis/tracking/TelescopeDetector.hpp new file mode 100644 index 0000000..19bb2b7 --- /dev/null +++ b/source/tdis/tracking/TelescopeDetector.hpp @@ -0,0 +1,64 @@ +// This file is part of the Acts project. +// +// Copyright (C) 2020 CERN for the benefit of the Acts project +// +// This Source Code Form is subject to the terms of the Mozilla Public +// License, v. 2.0. If a copy of the MPL was not distributed with this +// file, You can obtain one at http://mozilla.org/MPL/2.0/. + +#pragma once + +#include "Acts/Definitions/Units.hpp" +#include "Acts/Utilities/Logger.hpp" +#include "ActsExamples/Utilities/Options.hpp" + +#include +#include +#include +#include + +using namespace Acts::UnitLiterals; + +namespace Acts { +class TrackingGeometry; +class IMaterialDecorator; +} // namespace Acts + +namespace ActsExamples { +class IContextDecorator; +} // namespace ActsExamples + +namespace ActsExamples::Telescope { + +class TelescopeDetectorElement; +class TelescopeG4DetectorConstruction; + +struct TelescopeDetector { + using DetectorElement = ActsExamples::Telescope::TelescopeDetectorElement; + using DetectorElementPtr = std::shared_ptr; + using DetectorStore = std::vector; + + using ContextDecorators = + std::vector>; + using TrackingGeometryPtr = std::shared_ptr; + + struct Config { + std::vector positions{{0, 30, 60, 120, 150, 180}}; + std::vector stereos{{0, 0, 0, 0, 0, 0}}; + std::array offsets{{0, 0}}; + std::array bounds{{25, 100}}; + double thickness{80_um}; + int surfaceType{0}; + int binValue{2}; + }; + + Config config; + /// The store of the detector elements (lifetime: job) + DetectorStore detectorStore; + + std::pair finalize( + const Config& cfg, + const std::shared_ptr& mdecorator); +}; + +} // namespace ActsExamples::Telescope diff --git a/source/tdis/tracking/TelescopeDetectorElement.cpp b/source/tdis/tracking/TelescopeDetectorElement.cpp new file mode 100644 index 0000000..c761d43 --- /dev/null +++ b/source/tdis/tracking/TelescopeDetectorElement.cpp @@ -0,0 +1,42 @@ +// This file is part of the Acts project. +// +// Copyright (C) 2020-2021 CERN for the benefit of the Acts project +// +// This Source Code Form is subject to the terms of the Mozilla Public +// License, v. 2.0. If a copy of the MPL was not distributed with this +// file, You can obtain one at http://mozilla.org/MPL/2.0/. + +#include "TelescopeDetectorElement.hpp" + +#include "Acts/Surfaces/DiscSurface.hpp" +#include "Acts/Surfaces/PlaneSurface.hpp" + +ActsExamples::Telescope::TelescopeDetectorElement::TelescopeDetectorElement( + std::shared_ptr transform, + std::shared_ptr pBounds, double thickness, + std::shared_ptr material) + : Acts::DetectorElementBase(), + m_elementTransform(std::move(transform)), + m_elementSurface( + Acts::Surface::makeShared(pBounds, *this)), + m_elementThickness(thickness), + m_elementPlanarBounds(std::move(pBounds)), + m_elementDiscBounds(nullptr) { + auto mutableSurface = + std::const_pointer_cast(m_elementSurface); + mutableSurface->assignSurfaceMaterial(std::move(material)); +} + +ActsExamples::Telescope::TelescopeDetectorElement::TelescopeDetectorElement( + std::shared_ptr transform, + std::shared_ptr dBounds, double thickness, + std::shared_ptr material) + : Acts::DetectorElementBase(), + m_elementTransform(std::move(transform)), + m_elementSurface( + Acts::Surface::makeShared(dBounds, *this)), + m_elementThickness(thickness), + m_elementPlanarBounds(nullptr), + m_elementDiscBounds(std::move(dBounds)) { + m_elementSurface->assignSurfaceMaterial(std::move(material)); +} diff --git a/source/tdis/tracking/TelescopeDetectorElement.hpp b/source/tdis/tracking/TelescopeDetectorElement.hpp new file mode 100644 index 0000000..e56bbf9 --- /dev/null +++ b/source/tdis/tracking/TelescopeDetectorElement.hpp @@ -0,0 +1,162 @@ +// This file is part of the Acts project. +// +// Copyright (C) 2020-2021 CERN for the benefit of the Acts project +// +// This Source Code Form is subject to the terms of the Mozilla Public +// License, v. 2.0. If a copy of the MPL was not distributed with this +// file, You can obtain one at http://mozilla.org/MPL/2.0/. + +#pragma once + +#include "Acts/Definitions/Algebra.hpp" +#include "Acts/Geometry/DetectorElementBase.hpp" +#include "Acts/Geometry/GeometryContext.hpp" +#include "Acts/Surfaces/Surface.hpp" + +#include +#include +#include + +namespace Acts { +class Surface; +class PlanarBounds; +class DiscBounds; +class ISurfaceMaterial; +} // namespace Acts + +namespace ActsExamples::Telescope { + +/// @class TelescopeDetectorElement +/// +/// This is a lightweight type of detector element, +/// it simply implements the base class. +class TelescopeDetectorElement : public Acts::DetectorElementBase { + public: + /// @class ContextType + /// convention: nested to the Detector element + struct ContextType { + /// The current interval of validity + unsigned int iov = 0; + }; + + /// Constructor for single sided detector element + /// - bound to a Plane Surface + /// + /// @param transform is the transform that element the layer in 3D frame + /// @param pBounds is the planar bounds for the planar detector element + /// @param thickness is the module thickness + /// @param material is the (optional) Surface material associated to it + TelescopeDetectorElement( + std::shared_ptr transform, + std::shared_ptr pBounds, double thickness, + std::shared_ptr material = nullptr); + + /// Constructor for single sided detector element + /// - bound to a Disc Surface + /// + /// @param transform is the transform that element the layer in 3D frame + /// @param dBounds is the planar bounds for the disc like detector element + /// @param thickness is the module thickness + /// @param material is the (optional) Surface material associated to it + TelescopeDetectorElement( + std::shared_ptr transform, + std::shared_ptr dBounds, double thickness, + std::shared_ptr material = nullptr); + + /// Destructor + ~TelescopeDetectorElement() override = default; + + /// Return surface associated with this detector element + const Acts::Surface& surface() const final; + + /// Non-const access to the surface associated with this detector element + Acts::Surface& surface() final; + + /// The maximal thickness of the detector element wrt normal axis + double thickness() const final; + + /// Return local to global transform associated with this identifier + /// + /// @param gctx The current geometry context object, e.g. alignment + /// + /// @note this is called from the surface().transform() in the PROXY mode + const Acts::Transform3& transform( + const Acts::GeometryContext& gctx) const final; + + /// Return the nominal local to global transform + /// + /// @note the geometry context will hereby be ignored + const Acts::Transform3& nominalTransform( + const Acts::GeometryContext& gctx) const; + + /// Return local to global transform associated with this identifier + /// + /// @param alignedTransform is a new transform + /// @param iov is the batch for which it is meant + void addAlignedTransform(std::unique_ptr alignedTransform, + unsigned int iov); + + /// Return the set of alignment transforms in flight + const std::vector>& alignedTransforms() + const; + + private: + /// the transform for positioning in 3D space + std::shared_ptr m_elementTransform = nullptr; + // the aligned transforms + std::vector> m_alignedTransforms = {}; + /// the surface represented by it + std::shared_ptr m_elementSurface = nullptr; + /// the element thickness + double m_elementThickness = 0.; + /// the planar bounds + std::shared_ptr m_elementPlanarBounds = nullptr; + /// the disc bounds + std::shared_ptr m_elementDiscBounds = nullptr; +}; + +inline const Acts::Surface& TelescopeDetectorElement::surface() const { + return *m_elementSurface; +} + +inline Acts::Surface& TelescopeDetectorElement::surface() { + return *m_elementSurface; +} + +inline double TelescopeDetectorElement::thickness() const { + return m_elementThickness; +} + +inline const Acts::Transform3& TelescopeDetectorElement::transform( + const Acts::GeometryContext& gctx) const { + // Check if a different transform than the nominal exists + if (!m_alignedTransforms.empty()) { + // cast into the right context object + auto alignContext = gctx.get(); + return (*m_alignedTransforms[alignContext.iov].get()); + } + // Return the standard transform if not found + return nominalTransform(gctx); +} + +inline const Acts::Transform3& TelescopeDetectorElement::nominalTransform( + const Acts::GeometryContext& /*gctx*/) const { + return *m_elementTransform; +} + +inline void TelescopeDetectorElement::addAlignedTransform( + std::unique_ptr alignedTransform, unsigned int iov) { + // most standard case, it's just a new one: + auto cios = m_alignedTransforms.size(); + for (unsigned int ic = cios; ic <= iov; ++ic) { + m_alignedTransforms.push_back(nullptr); + } + m_alignedTransforms[iov] = std::move(alignedTransform); +} + +inline const std::vector>& +TelescopeDetectorElement::alignedTransforms() const { + return m_alignedTransforms; +} + +} // namespace ActsExamples::Telescope From ca1ea9818dae85c9f1da010530c8028b00b063eb Mon Sep 17 00:00:00 2001 From: Dmitry Romanov Date: Fri, 15 Nov 2024 11:31:22 -0500 Subject: [PATCH 4/8] Geometry export to OBJ and PLY --- source/tdis/tracking/ActsGeometryService.cc | 35 ++++++---- source/tdis/tracking/ActsGeometryService.h | 66 ++++++++++++------- .../tdis/tracking/BuildTelescopeDetector.cpp | 4 +- 3 files changed, 70 insertions(+), 35 deletions(-) diff --git a/source/tdis/tracking/ActsGeometryService.cc b/source/tdis/tracking/ActsGeometryService.cc index 4f21c57..6ce18a7 100644 --- a/source/tdis/tracking/ActsGeometryService.cc +++ b/source/tdis/tracking/ActsGeometryService.cc @@ -16,6 +16,7 @@ #include #include #include +#include #include #include #include @@ -357,7 +358,7 @@ void tdis::tracking::ActsGeometryService::Init() { // We just copy it, considering the disks are the same bounds[0] = r_min * Acts::UnitConstants::cm; bounds[1] = r_max * Acts::UnitConstants::cm; - thickness = dz * Acts::UnitConstants::cm * 0.1; + thickness = dz * Acts::UnitConstants::cm * 0.001; m_log->info(fmt::format(" Position: ({:.4f}, {:.4f}, {:>8.4f}) cm, rmin: {:.3f} cm, rmax: {:.3f} cm, dz: {:.4f} cm, Node: '{}', Volume: '{}'", x, y, z, r_min, r_max, dz, node->GetName(), volume->GetName())); @@ -366,14 +367,14 @@ void tdis::tracking::ActsGeometryService::Init() { stereos.push_back(0); } - Acts::ObjVisualization3D objVis; + m_init_log->info("Building ACTS Geometry"); /// Return the telescope detector - std::shared_ptr gGeometry = + gGeometry = ActsExamples::Telescope::buildDetector( nominalContext, // gctx is the detector element dependent geometry context detectorStore, // detectorStore is the store for the detector element @@ -385,20 +386,32 @@ void tdis::tracking::ActsGeometryService::Init() { ActsExamples::Telescope::TelescopeSurfaceType::Disc, // surfaceType is the detector surface type Acts::BinningValue::binZ); - m_log->info("ACTS Detector Elements"); - for(auto &element: detectorStore) { - Acts::GeometryView3D::drawSurface(objVis, element->surface(), m_geometry_context, element->transform(m_geometry_context)); - } + // Visualize ACTS geometry + const Acts::TrackingVolume& tgVolume = *(gGeometry->highestTrackingVolume()); + + // OBJ visualization export + auto obj_file = m_obj_output_file(); if(!m_obj_output_file().empty()) { - m_log->info("ACTS Exporting geometry as view to: {}", m_obj_output_file()); - objVis.write(m_obj_output_file()); + m_log->info("ACTS exporting to OBJ: {}", m_obj_output_file()); + Acts::ObjVisualization3D vis_helper; + DrawTrackingGeometry(vis_helper, tgVolume, m_geometry_context); + vis_helper.write(m_obj_output_file()); } else { - m_log->info("ACTS geometry filename defined by: '{}' is empty. NOT exporting geometry to .obj or whatever ", m_obj_output_file.m_name); + m_log->info("ACTS OBJ Export. Flag '{}' is empty. NOT exporting.", m_obj_output_file.m_name); } - exit(0); + // PLY visualization export + if(!m_ply_output_file().empty()) { + m_log->info("ACTS exporting to PLY: {}", m_ply_output_file()); + Acts::PlyVisualization3D vis_helper; + DrawTrackingGeometry(vis_helper, tgVolume, m_geometry_context); + vis_helper.write(m_ply_output_file()); + } else { + m_log->info("ACTS PLY Export. Flag '{}' is empty. NOT exporting.", m_ply_output_file.m_name); + } + exit(0); // Set ticker back m_app->SetTicker(was_ticker_enabled); diff --git a/source/tdis/tracking/ActsGeometryService.h b/source/tdis/tracking/ActsGeometryService.h index d646766..3f2284a 100644 --- a/source/tdis/tracking/ActsGeometryService.h +++ b/source/tdis/tracking/ActsGeometryService.h @@ -9,46 +9,41 @@ #include #include -#include "TelescopeDetectorElement.hpp" - +#include #include #include #include "ActsGeometryProvider.h" +#include "TelescopeDetectorElement.hpp" #include "services/LogService.hpp" namespace tdis::tracking { - class ActsGeometryService : public JService - { - public: + class ActsGeometryService : public JService { + public: explicit ActsGeometryService() : JService() {} ~ActsGeometryService() override = default; - void Init() override; TGeoManager* GetGeoManager() const { return m_tgeo_manager; } - TGeoVolume * GetTopVolume() const { return m_tgeo_manager->GetTopVolume(); } - - TGeoNode * GetTopNode() const { return m_tgeo_manager->GetTopNode(); } + TGeoVolume* GetTopVolume() const { return m_tgeo_manager->GetTopVolume(); } + TGeoNode* GetTopNode() const { return m_tgeo_manager->GetTopNode(); } - Acts::GeometryContext& GetActsGeometryContext() { - return m_geometry_context; - } + Acts::GeometryContext& GetActsGeometryContext() { return m_geometry_context; } + template void DrawTrackingGeometry(TVis& vis_helper, + const Acts::TrackingVolume& tVolume, + const Acts::GeometryContext& gctx); - private: - - - Parameter m_tgeo_file {this, "acts:geometry", "g4sbs_mtpc.root", "TGeo filename with geometry for ACTS"}; - Parameter m_material_map_file {this, "acts:material_map", "", "JSON/CBOR material map file path"}; - Parameter m_obj_output_file {this, "acts:output_obj", "", "Output file name to dump ACTS converted geometry as OBJ"}; - Parameter m_ply_output_file {this, "acts:output_ply", "", "Output file name to dump ACTS converted geometry as PLY"}; - - Service m_service_log {this}; + private: + Parameter m_tgeo_file{this, "acts:geometry", "g4sbs_mtpc.root","TGeo filename with geometry for ACTS"}; + Parameter m_material_map_file{this, "acts:material_map", "", "JSON/CBOR material map file path"}; + Parameter m_obj_output_file{this, "acts:output_obj", "", "Output file name to dump ACTS converted geometry as OBJ"}; + Parameter m_ply_output_file{this, "acts:output_ply", "", "Output file name to dump ACTS converted geometry as PLY"}; + Service m_service_log{this}; // General acts log std::shared_ptr m_log; @@ -65,7 +60,34 @@ namespace tdis::tracking { ActsExamples::Telescope::TelescopeDetectorElement::ContextType nominalContext; - std::vector> detectorStore; + std::vector> + detectorStore; + std::shared_ptr gGeometry; }; + + + template + void ActsGeometryService::DrawTrackingGeometry(TVis& vis_helper, const Acts::TrackingVolume& tVolume, const Acts::GeometryContext& gctx) { + bool triangulate = false; + + Acts::ViewConfig viewSensitive{.color = {0, 180, 240}, .quarterSegments = 72, .triangulate = triangulate}; + Acts::ViewConfig viewPassive{.color = {240, 280, 0}, .quarterSegments = 72, .triangulate = triangulate}; + Acts::ViewConfig viewVolume{.color = {220, 220, 0}, .quarterSegments = 72, .triangulate = triangulate}; + Acts::ViewConfig viewContainer{.color = {220, 220, 0}, .quarterSegments = 72, .triangulate = triangulate}; + Acts::ViewConfig viewGrid{.color = {220, 0, 0}, .offset = 3., .quarterSegments = 8, .triangulate = triangulate}; + + Acts::GeometryView3D::drawTrackingVolume( + vis_helper, // Visualization helper (templated) + tVolume, // Tracking volume to be drawn + gctx, // Geometry context + viewContainer, // Container volume view configuration + viewVolume, // Navigation level volume view configuration + viewPassive, // Passive surfaces view configuration + viewSensitive, // Sensitive surfaces view configuration + viewGrid, // Grid display view configuration + true, // Write or not + "tdis" // Optional tag + ); + } } // namespace tdis::tracking \ No newline at end of file diff --git a/source/tdis/tracking/BuildTelescopeDetector.cpp b/source/tdis/tracking/BuildTelescopeDetector.cpp index c237127..3f7835c 100644 --- a/source/tdis/tracking/BuildTelescopeDetector.cpp +++ b/source/tdis/tracking/BuildTelescopeDetector.cpp @@ -90,9 +90,9 @@ ActsExamples::Telescope::buildDetector(const typename ActsExamples::Telescope::T std::unique_ptr surArray(new Acts::SurfaceArray(surface)); // Construct the layer if (surfaceType == TelescopeSurfaceType::Plane) { - layers[i] = Acts::PlaneLayer::create(trafo, pBounds, std::move(surArray), 1._mm); + layers[i] = Acts::PlaneLayer::create(trafo, pBounds, std::move(surArray), thickness); } else { - layers[i] = Acts::DiscLayer::create(trafo, rBounds, std::move(surArray), 1._mm); + layers[i] = Acts::DiscLayer::create(trafo, rBounds, std::move(surArray), thickness); } // Associate the layer to the surface auto mutableSurface = const_cast(surface.get()); From 4428f066c29bfb2b18f218f2a7ba5760ae113c25 Mon Sep 17 00:00:00 2001 From: Dmitry Romanov Date: Fri, 15 Nov 2024 23:30:39 -0500 Subject: [PATCH 5/8] Reconstruction tiests --- .github/workflows/linux-build-test.yml | 2 +- README.md | 32 +-- source/tdis/CMakeLists.txt | 2 - source/tdis/PadGeometryHelper.hpp | 25 ++- source/tdis/io/DigitizedDataEventSource.hpp | 15 +- source/tdis/io/PodioWriteProcessor.hpp | 5 +- source/tdis/layout.yaml | 18 +- source/tdis/tdis_main.cpp | 2 +- source/tdis/tracking/ActsGeometryProvider.cc | 194 ------------------ source/tdis/tracking/ActsGeometryProvider.h | 142 ------------- source/tdis/tracking/ActsGeometryService.cc | 27 +-- source/tdis/tracking/ActsGeometryService.h | 9 +- .../tdis/tracking/ReconstructedHitFactory.h | 140 ++++++------- 13 files changed, 131 insertions(+), 482 deletions(-) delete mode 100644 source/tdis/tracking/ActsGeometryProvider.cc delete mode 100644 source/tdis/tracking/ActsGeometryProvider.h diff --git a/.github/workflows/linux-build-test.yml b/.github/workflows/linux-build-test.yml index 18cd491..c26f6e0 100644 --- a/.github/workflows/linux-build-test.yml +++ b/.github/workflows/linux-build-test.yml @@ -28,7 +28,7 @@ jobs: key: ${{ github.workflow }}-cpm-modules-${{ hashFiles('**/CMakeLists.txt', '**/*.cmake') }} - name: configure - run: cmake -Stest -Bbuild -DENABLE_TEST_COVERAGE=1 -DCMAKE_BUILD_TYPE=Debug + run: cmake -Bbuild -DCMAKE_BUILD_TYPE=Debug - name: build run: cmake --build build -j4 diff --git a/README.md b/README.md index 8d3a176..84e7725 100644 --- a/README.md +++ b/README.md @@ -12,8 +12,9 @@ tdis /mnt/data/g4sbsout_EPCEvents_200000.txt ``` +## -# Geometry IDs +## Geometry IDs The geometry description is gives us that each pad is identified by `plane`, `ring`, `pad`: @@ -21,29 +22,6 @@ The geometry description is gives us that each pad is identified by `plane`, `ri - `ring` The pads are arranged into 21 concentric rings, The rings are numbered from 0 to 20 sequentially, with 0 being the innermost ring. - `pad` Each ring has 122 pads, 0 is at or closest to phi=0 and numbering is clockwise -This is convenient to encode in a singe uint32_t ID. We assign 8 bits each to plane and ring, and 16 bits to pad. -The encodePadID function packs these three values into a single uint32_t using bitwise shifts and OR operations. - -```python -import numpy as np - -def encode_pad_id(plane, ring, pad): - """Encode plane, ring, pad into an uint32_t id. - Works with both scalar values and numpy arrays.""" - plane = np.asarray(plane, dtype=np.uint32) - ring = np.asarray(ring, dtype=np.uint32) - pad = np.asarray(pad, dtype=np.uint32) - - pad_id = (plane << 24) | (ring << 16) | pad - return pad_id - -def decode_pad_id(pad_id): - """Decode pad_id into plane, ring, pad. - Works with both scalar values and numpy arrays.""" - - pad_id = np.asarray(pad_id, dtype=np.uint32) - plane = (pad_id >> 24) & 0xFF # Upper 8 bits - ring = (pad_id >> 16) & 0xFF # Next 8 bits - pad = pad_id & 0xFFFF # Lower 16 bits - return plane, ring, pad -``` \ No newline at end of file +This is convenient to encode in a singe uint32_t ID. + +`1 000 000 * plane + 1 000 * ring + pad` \ No newline at end of file diff --git a/source/tdis/CMakeLists.txt b/source/tdis/CMakeLists.txt index 5faa935..f6d3e66 100644 --- a/source/tdis/CMakeLists.txt +++ b/source/tdis/CMakeLists.txt @@ -40,8 +40,6 @@ add_executable(tdis PadGeometryHelper.hpp io/DigitizedDataEventSource.hpp io/PodioWriteProcessor.hpp - tracking/ActsGeometryProvider.cc - tracking/ActsGeometryProvider.h tracking/ActsGeometryService.cc tracking/ActsGeometryService.h tracking/ReconstructedHitFactory.h diff --git a/source/tdis/PadGeometryHelper.hpp b/source/tdis/PadGeometryHelper.hpp index 8d24484..ad7261d 100644 --- a/source/tdis/PadGeometryHelper.hpp +++ b/source/tdis/PadGeometryHelper.hpp @@ -5,25 +5,37 @@ #include // For std::pair // Constants -constexpr int num_rings = 21; +constexpr int num_rings = 21 * Acts::UnitConstants::cm; constexpr int num_pads_per_ring = 122; -constexpr double min_radius = 5.0; // Minimum radius in cm -constexpr double max_radius = 15.0; // Maximum radius in cm +constexpr double min_radius = 5.0 * Acts::UnitConstants::cm; // Minimum radius in cm +constexpr double max_radius = 15.0 * Acts::UnitConstants::cm; // Maximum radius in cm constexpr double total_radius = max_radius - min_radius; // Total radial width (10 cm) constexpr double ring_width = total_radius / num_rings; // Radial width of each ring constexpr double delta_theta = 2 * M_PI / num_pads_per_ring; // Angular width of each pad in radians namespace tdis { - inline std::tuple getPadCenter(int ring, int pad) { + + inline double getPadHight() { + /** gets pad height which is the distance between rings*/ + return (max_radius - min_radius) / num_pads_per_ring; + } + + inline double getPadApproxWidth(const int ring) { + /** gets pad approximate width calculated from curvature, which is not exact the width */ + const double r = min_radius + (ring + 0.5) * getPadHight(); + return r * 2 * M_PI / num_pads_per_ring; + } + + inline std::tuple getPadCenter(const int ring, const int pad) { /* Compute the X and Y coordinates of the center of a pad given its ring and pad indices. Parameters: - ring (int): The ring index (0 is the innermost ring). - - pad (int): The pad index (0 is at or closest to φ=0°, numbering is clockwise). + - pad (int): The pad index (0 is at or closest to φ=0, numbering is clockwise). Returns: - - std::pair: X and Y coordinates of the pad center in cm. + - std::pair: X and Y coordinates of the pad center in Acts::UnitConstants. */ // Validate inputs @@ -57,4 +69,5 @@ namespace tdis { return { x, y }; } + } \ No newline at end of file diff --git a/source/tdis/io/DigitizedDataEventSource.hpp b/source/tdis/io/DigitizedDataEventSource.hpp index d810a7e..933956a 100644 --- a/source/tdis/io/DigitizedDataEventSource.hpp +++ b/source/tdis/io/DigitizedDataEventSource.hpp @@ -41,10 +41,9 @@ #include #include +#include #include -#include #include -#include #include #include #include @@ -345,13 +344,21 @@ namespace tdis::io { podioTrack.vertexZ(track.vertexZ); podioTrack.momentum(track.momentum); for(auto& hit: track.hits) { + auto podioHit = podioHits.create(); - podioHit.time( hit.time ); + podioHit.time( hit.time * Acts::UnitConstants::ns ); podioHit.adc( hit.adc ); podioHit.ring( hit.ring ); podioHit.pad( hit.pad ); podioHit.plane( hit.plane ); - podioHit.zToGem( hit.zToGem); + podioHit.zToGem( hit.zToGem * Acts::UnitConstants::m); + + edm4hep::Vector3f true_pos = edm4hep::Vector3f{ + static_cast(hit.true_x * Acts::UnitConstants::m), + static_cast(hit.true_y * Acts::UnitConstants::m), + static_cast(hit.true_z * Acts::UnitConstants::m) + }; + podioHit.truePosition(true_pos); podioTrack.addhits(podioHit); } diff --git a/source/tdis/io/PodioWriteProcessor.hpp b/source/tdis/io/PodioWriteProcessor.hpp index c060335..690a4be 100644 --- a/source/tdis/io/PodioWriteProcessor.hpp +++ b/source/tdis/io/PodioWriteProcessor.hpp @@ -75,7 +75,10 @@ class PodioWriteProcessor : public JEventProcessor { // Truth record "DigitizedMtpcMcTrack", - "DigitizedMtpcMcHit" + "DigitizedMtpcMcHit", + + // Digitized hits + "TrackerHit" }; PodioWriteProcessor(JApplication * app); diff --git a/source/tdis/layout.yaml b/source/tdis/layout.yaml index 4c2c491..6bb2368 100644 --- a/source/tdis/layout.yaml +++ b/source/tdis/layout.yaml @@ -283,15 +283,13 @@ datatypes : Description: "TDIS MTPC Digitized track" Author: "Dmitry Romanov" Members: - - double time // - Time of arrival at Pad (ns) - - double adc // - Amplitude (ADC bin of sampa) - - int ring // - Ring (id of rin, 0 is inner most). - - int pad // - Pad (id of pad, 0 is at or closest to phi=0 and numbering is clockwise). - - int plane // - Plane(id of z plane from 0 upstream to 9 downstream) - - double zToGem // - ZtoGEM (m) - - double true_x // - True hit x info (quiet_NaN() if not provided) - - double true_y // - True hit y info (quiet_NaN() if not provided) - - double true_z // - True hit z info (quiet_NaN() if not provided) + - double time // - Time of arrival at Pad (ns) + - double adc // - Amplitude (ADC bin of sampa) + - int ring // - Ring (id of rin, 0 is inner most). + - int pad // - Pad (id of pad, 0 is at or closest to phi=0 and numbering is clockwise). + - int plane // - Plane(id of z plane from 0 upstream to 9 downstream) + - double zToGem // - ZtoGEM (m) + - edm4hep::Vector3f truePosition // - True hit x,y,z info (quiet_NaN() if not provided) edm4eic::TrackerHit: Description: "Tracker hit (reconstructed from Raw)" @@ -299,7 +297,7 @@ datatypes : Members: - uint64_t cellID // The detector specific (geometrical) cell id. - edm4hep::Vector3f position // Hit (cell) position [mm] - - edm4eic::CovDiag3f positionError // Covariance Matrix + - edm4eic::CovDiag3f positionError // Covariance Matrix - float time // Hit time [ns] - float timeError // Error on the time - float edep // Energy deposit in this hit [GeV] diff --git a/source/tdis/tdis_main.cpp b/source/tdis/tdis_main.cpp index 8954c66..601e599 100644 --- a/source/tdis/tdis_main.cpp +++ b/source/tdis/tdis_main.cpp @@ -116,7 +116,7 @@ int main(int argc, char* argv[]) { app.Add(new JEventSourceGeneratorT); app.Add(new tdis::io::PodioWriteProcessor(&app)); - app.Add(new JFactoryGeneratorT); + app.Add(new JFactoryGeneratorT("TrackerHit")); // app.Add(new JEventProcessorPodio); // app.Add(new JFactoryGeneratorT()); // app.Add(new JFactoryGeneratorT()); diff --git a/source/tdis/tracking/ActsGeometryProvider.cc b/source/tdis/tracking/ActsGeometryProvider.cc deleted file mode 100644 index 15d3e48..0000000 --- a/source/tdis/tracking/ActsGeometryProvider.cc +++ /dev/null @@ -1,194 +0,0 @@ -// SPDX-License-Identifier: LGPL-3.0-or-later -// Copyright (C) 2022 - 2024 Whitney Armstrong, Wouter Deconinck, Dmitry Romanov - -#include -#include -#include -#include -#include -#include - -#include -#include -#include -#include -#include -#include -#include -#include - -#include -#include -#include -#include -#include -#include -#include -#include - -#include "ActsGeometryProvider.h" -#include "extensions/spdlog/SpdlogToActs.h" - -// Formatter for Eigen matrices -#if FMT_VERSION >= 90000 -#include - -template -struct fmt::formatter< - T, - std::enable_if_t< - std::is_base_of_v, T>, - char - > -> : fmt::ostream_formatter {}; -#endif // FMT_VERSION >= 90000 - -void ActsGeometryProvider::initialize(const std::string& geometry_file, - const std::string& material_file, - std::shared_ptr log, - std::shared_ptr init_log) { - // LOGGING - m_log = log; - m_init_log = init_log; - - m_init_log->debug("ActsGeometryProvider initialization"); - - m_init_log->debug("Set TGeoManager and acts_init_log_level log levels"); - // Turn off TGeo printouts if appropriate for the msg level - if (m_log->level() >= (int) spdlog::level::info) { - TGeoManager::SetVerboseLevel(0); - } - - // Set ACTS logging level - auto acts_init_log_level = eicrecon::SpdlogToActsLevel(m_init_log->level()); - - // Load ACTS materials maps - std::shared_ptr materialDeco{nullptr}; - if (!material_file.empty()) { - m_init_log->info("loading materials map from file: '{}'", material_file); - // Set up the converter first - Acts::MaterialMapJsonConverter::Config jsonGeoConvConfig; - // Set up the json-based decorator - materialDeco = std::make_shared(jsonGeoConvConfig, material_file, acts_init_log_level); - } - - // // Geometry identifier hook to write detector ID to extra field - // TODO decide if we need this - // class ConvertDD4hepDetectorGeometryIdentifierHook: public Acts::GeometryIdentifierHook { - // Acts::GeometryIdentifier decorateIdentifier(Acts::GeometryIdentifier identifier, const Acts::Surface& surface) const { - // const auto* dd4hep_det_element = dynamic_cast(surface.associatedDetectorElement()); - // if (dd4hep_det_element == nullptr) { - // return identifier; - // } - // // set 8-bit extra field to 8-bit DD4hep detector ID - // return identifier.setExtra(0xff & dd4hep_det_element->identifier()); - // }; - // }; - // auto geometryIdHook = std::make_shared(); - - // Convert DD4hep geometry to ACTS - m_init_log->info("Converting DD4Hep geometry to ACTS..."); - auto logger = eicrecon::getSpdlogLogger("CONV", m_log); - // Acts::BinningType bTypePhi = Acts::equidistant; - // Acts::BinningType bTypeR = Acts::equidistant; - // Acts::BinningType bTypeZ = Acts::equidistant; - // double layerEnvelopeR = Acts::UnitConstants::mm; - // double layerEnvelopeZ = Acts::UnitConstants::mm; - // double defaultLayerThickness = Acts::UnitConstants::fm; - // using Acts::sortDetElementsByID; - // - // try { - // m_trackingGeo = Acts::convertDD4hepDetector( - // m_dd4hepDetector->world(), - // *logger, - // bTypePhi, - // bTypeR, - // bTypeZ, - // layerEnvelopeR, - // layerEnvelopeZ, - // defaultLayerThickness, - // sortDetElementsByID, - // m_trackingGeoCtx, - // materialDeco, - // geometryIdHook); - // } - // catch(std::exception &ex) { - // m_init_log->error("Error during DD4Hep -> ACTS geometry conversion: {}", ex.what()); - // m_init_log->info ("Set parameter acts::InitLogLevel=trace to see conversion info and possibly identify failing geometry"); - // throw JException(ex.what()); - // } - // - // m_init_log->info("DD4Hep geometry converted!"); - - // Visit surfaces - // TODO 2 Revisit this surface checking - // m_init_log->info("Checking surfaces..."); - // if (m_trackingGeo) { - // // Write tracking geometry to collection of obj or ply files - // const Acts::TrackingVolume* world = m_trackingGeo->highestTrackingVolume(); - // if (m_objWriteIt) { - // m_init_log->info("Writing obj files to {}...", m_outputDir); - // Acts::ObjVisualization3D objVis; - // Acts::GeometryView3D::drawTrackingVolume( - // objVis, *world, m_trackingGeoCtx, - // m_containerView, m_volumeView, m_passiveView, m_sensitiveView, m_gridView, - // m_objWriteIt, m_outputTag, m_outputDir - // ); - // } - // if (m_plyWriteIt) { - // m_init_log->info("Writing ply files to {}...", m_outputDir); - // Acts::PlyVisualization3D plyVis; - // Acts::GeometryView3D::drawTrackingVolume( - // plyVis, *world, m_trackingGeoCtx, - // m_containerView, m_volumeView, m_passiveView, m_sensitiveView, m_gridView, - // m_plyWriteIt, m_outputTag, m_outputDir - // ); - // } - // - // m_init_log->debug("visiting all the surfaces "); - // m_trackingGeo->visitSurfaces([this](const Acts::Surface *surface) { - // // for now we just require a valid surface - // if (surface == nullptr) { - // m_init_log->info("no surface??? "); - // return; - // } - // const auto *det_element = - // dynamic_cast(surface->associatedDetectorElement()); - // - // if (det_element == nullptr) { - // m_init_log->error("invalid det_element!!! det_element == nullptr "); - // return; - // } - // - // // more verbose output is lower enum value - // m_init_log->debug(" det_element->identifier() = {} ", det_element->identifier()); - // auto volman = m_dd4hepDetector->volumeManager(); - // auto *vol_ctx = volman.lookupContext(det_element->identifier()); - // auto vol_id = vol_ctx->identifier; - // - // if (m_init_log->level() <= spdlog::level::debug) { - // auto de = vol_ctx->element; - // m_init_log->debug(" de.path = {}", de.path()); - // m_init_log->debug(" de.placementPath = {}", de.placementPath()); - // } - // - // this->m_surfaces.insert_or_assign(vol_id, surface); - // }); - // } - // else { - // m_init_log->error("m_trackingGeo==null why am I still alive???"); - // } - - // Load ACTS magnetic field - // TODO we don't need this - // m_init_log->info("Loading magnetic field..."); - // m_magneticField = std::make_shared(m_dd4hepDetector); - // Acts::MagneticFieldContext m_fieldctx{eicrecon::BField::BFieldVariant(m_magneticField)}; - // auto bCache = m_magneticField->makeCache(m_fieldctx); - // for (int z: {0, 500, 1000, 1500, 2000, 3000, 4000}) { - // auto b = m_magneticField->getField({0.0, 0.0, double(z)}, bCache).value(); - // m_init_log->debug("B(z = {:>5} [mm]) = {} T", z, b.transpose() / Acts::UnitConstants::T); - // } - - m_init_log->info("ActsGeometryProvider initialization complete"); -} diff --git a/source/tdis/tracking/ActsGeometryProvider.h b/source/tdis/tracking/ActsGeometryProvider.h deleted file mode 100644 index 6d4b32b..0000000 --- a/source/tdis/tracking/ActsGeometryProvider.h +++ /dev/null @@ -1,142 +0,0 @@ -// SPDX-License-Identifier: LGPL-3.0-or-later -// Copyright (C) 2015 - 2024 Julia Hrdinka, Whitney Armstrong, Wouter Deconinck, Dmitry Romanov - -#pragma once - -// ACTS -#include -#include -#include -#include -#include -#include -#include -#include -#include -#include -#include -#include -#include -#include -#include -#include - -#include "DD4hepBField.h" - -namespace dd4hep::rec { - class Surface; -} - -/** Draw the surfaces and save to obj file. - * This is useful for debugging the ACTS geometry. The obj file can - * be loaded into various tools, such as FreeCAD, for inspection. - */ -void draw_surfaces(std::shared_ptr trk_geo, std::shared_ptr init_log, const std::string &fname); - -class ActsGeometryProvider { -public: - ActsGeometryProvider() {} - using VolumeSurfaceMap = std::unordered_map; - - virtual void initialize(const std::string& geometry_file, - const std::string& material_file, - std::shared_ptr log, - std::shared_ptr init_log) final; - - - /** Gets the ACTS tracking geometry. - */ - std::shared_ptr trackingGeometry() const { return m_trackingGeo;} - - std::shared_ptr getFieldProvider() const { return m_magneticField; } - - - const VolumeSurfaceMap &surfaceMap() const { return m_surfaces; } - - - std::map getDD4hepSurfaceMap() const { return m_surfaceMap; } - - const Acts::GeometryContext& getActsGeometryContext() const {return m_trackingGeoCtx;} - - /// ACTS general logger that is used for running ACTS - std::shared_ptr getActsRelatedLogger() const { return m_log; } - - /// Logger that is used for geometry initialization - /// By default its level the same as ACTS general logger (m_log) - /// But it might be customized to solely printout geometry information - std::shared_ptr getActsInitRelatedLogger() const { return m_init_log; } - -private: - - - /// DD4hep surface map - std::map m_surfaceMap; - - /// ACTS Logging Level - Acts::Logging::Level acts_log_level = Acts::Logging::INFO; - - /// ACTS Tracking Geometry Context - Acts::GeometryContext m_trackingGeoCtx; - - /// ACTS Tracking Geometry - std::shared_ptr m_trackingGeo{nullptr}; - - /// ACTS surface lookup container for hit surfaces that generate smeared hits - VolumeSurfaceMap m_surfaces; - - /// Acts magnetic field - std::shared_ptr m_magneticField = nullptr; - - /// ACTS general logger that is used for running ACTS - std::shared_ptr m_log; - - /// Logger that is used for geometry initialization - /// By default its level the same as ACTS general logger (m_log) - /// But it might be customized to solely printout geometry information - std::shared_ptr m_init_log; - - /// Configuration for obj export - Acts::ViewConfig m_containerView{true, {220, 220, 220}}; - Acts::ViewConfig m_volumeView{true, {220, 220, 0}}; - Acts::ViewConfig m_sensitiveView{true, {0, 180, 240}}; - Acts::ViewConfig m_passiveView{true, {240, 280, 0}}; - Acts::ViewConfig m_gridView{true, {220, 0, 0}}; - bool m_objWriteIt{false}; - bool m_plyWriteIt{false}; - std::string m_outputTag{""}; - std::string m_outputDir{""}; - -public: - void setObjWriteIt(bool writeit) { m_objWriteIt = writeit; } - bool getObjWriteIt() const { return m_objWriteIt; } - void setPlyWriteIt(bool writeit) { m_plyWriteIt = writeit; } - bool getPlyWriteIt() const { return m_plyWriteIt; } - - void setOutputTag(std::string tag) { m_outputTag = tag; } - std::string getOutputTag() const { return m_outputTag; } - void setOutputDir(std::string dir) { m_outputDir = dir; } - std::string getOutputDir() const { return m_outputDir; } - - void setContainerView(std::array view) { - m_containerView.color = Acts::Color{view}; - } - const Acts::ViewConfig& getContainerView() const { return m_containerView; } - void setVolumeView(std::array view) { - m_volumeView.color = Acts::Color{view}; - } - const Acts::ViewConfig& getVolumeView() const { return m_volumeView; } - void setSensitiveView(std::array view) { - m_sensitiveView.color = Acts::Color{view}; - } - const Acts::ViewConfig& getSensitiveView() const { return m_sensitiveView; } - void setPassiveView(std::array view) { - m_passiveView.color = Acts::Color{view}; - } - - const Acts::ViewConfig& getPassiveView() const { return m_passiveView; } - void setGridView(std::array view) { - m_gridView.color = Acts::Color{view}; - } - const Acts::ViewConfig& getGridView() const { return m_gridView; } - -}; diff --git a/source/tdis/tracking/ActsGeometryService.cc b/source/tdis/tracking/ActsGeometryService.cc index 6ce18a7..8000952 100644 --- a/source/tdis/tracking/ActsGeometryService.cc +++ b/source/tdis/tracking/ActsGeometryService.cc @@ -23,7 +23,6 @@ #include #include -#include "ActsGeometryProvider.h" #include "BuildTelescopeDetector.hpp" #include "TelescopeDetector.hpp" #include "services/LogService.hpp" @@ -324,15 +323,11 @@ void tdis::tracking::ActsGeometryService::Init() { findNodesWithPrefix(m_tgeo_manager->GetTopNode(), "mTPCReadoutDisc", disk_nodes); - - - std::vector positions; + m_plane_positions.clear(); std::vector stereos; std::array offsets{{0, 0}}; std::array bounds{{25, 100}}; double thickness{0.1}; - int surfaceType{(int) ActsExamples::Telescope::TelescopeSurfaceType::Disc}; - m_init_log->info("Found readout disks: {}", disk_nodes.size()); for (auto node : disk_nodes) { @@ -350,7 +345,7 @@ void tdis::tracking::ActsGeometryService::Init() { continue; } - TGeoTube* tube = dynamic_cast(shape); + auto* tube = dynamic_cast(shape); double r_min = tube->GetRmin(); double r_max = tube->GetRmax(); double dz = tube->GetDz(); @@ -363,7 +358,7 @@ void tdis::tracking::ActsGeometryService::Init() { m_log->info(fmt::format(" Position: ({:.4f}, {:.4f}, {:>8.4f}) cm, rmin: {:.3f} cm, rmax: {:.3f} cm, dz: {:.4f} cm, Node: '{}', Volume: '{}'", x, y, z, r_min, r_max, dz, node->GetName(), volume->GetName())); - positions.push_back(z * Acts::UnitConstants::cm); + m_plane_positions.push_back(z * Acts::UnitConstants::cm); stereos.push_back(0); } @@ -376,13 +371,13 @@ void tdis::tracking::ActsGeometryService::Init() { /// Return the telescope detector gGeometry = ActsExamples::Telescope::buildDetector( - nominalContext, // gctx is the detector element dependent geometry context - detectorStore, // detectorStore is the store for the detector element - positions, // positions are the positions of different layers in the longitudinal direction - stereos, // stereoAngles are the stereo angles of different layers, which are rotation angles around the longitudinal (normal) direction - offsets, // is the offset (u, v) of the layers in the transverse plane - bounds, // bounds is the surface bound values, i.e. halfX and halfY if plane surface, and minR and maxR if disc surface - thickness, // thickness is the material thickness of each layer + nominalContext, // gctx is the detector element dependent geometry context + detectorStore, // detectorStore is the store for the detector element + m_plane_positions, // positions are the positions of different layers in the longitudinal direction + stereos, // stereoAngles are the stereo angles of different layers, which are rotation angles around the longitudinal (normal) direction + offsets, // is the offset (u, v) of the layers in the transverse plane + bounds, // bounds is the surface bound values, i.e. halfX and halfY if plane surface, and minR and maxR if disc surface + thickness, // thickness is the material thickness of each layer ActsExamples::Telescope::TelescopeSurfaceType::Disc, // surfaceType is the detector surface type Acts::BinningValue::binZ); @@ -411,8 +406,6 @@ void tdis::tracking::ActsGeometryService::Init() { m_log->info("ACTS PLY Export. Flag '{}' is empty. NOT exporting.", m_ply_output_file.m_name); } - exit(0); - // Set ticker back m_app->SetTicker(was_ticker_enabled); diff --git a/source/tdis/tracking/ActsGeometryService.h b/source/tdis/tracking/ActsGeometryService.h index 3f2284a..fe99e27 100644 --- a/source/tdis/tracking/ActsGeometryService.h +++ b/source/tdis/tracking/ActsGeometryService.h @@ -9,11 +9,11 @@ #include #include +#include #include #include #include -#include "ActsGeometryProvider.h" #include "TelescopeDetectorElement.hpp" #include "services/LogService.hpp" @@ -33,6 +33,8 @@ namespace tdis::tracking { Acts::GeometryContext& GetActsGeometryContext() { return m_geometry_context; } + const std::vector& GetPlanePositions() const { return m_plane_positions; } + template void DrawTrackingGeometry(TVis& vis_helper, const Acts::TrackingVolume& tVolume, const Acts::GeometryContext& gctx); @@ -54,7 +56,7 @@ namespace tdis::tracking { std::shared_ptr m_init_log; /// Root TGeo Manater for TGeo Geometry - TGeoManager* m_tgeo_manager; + TGeoManager* m_tgeo_manager = nullptr; Acts::GeometryContext m_geometry_context = Acts::GeometryContext(); @@ -64,6 +66,9 @@ namespace tdis::tracking { detectorStore; std::shared_ptr gGeometry; + + // Plane positions + std::vector m_plane_positions; }; diff --git a/source/tdis/tracking/ReconstructedHitFactory.h b/source/tdis/tracking/ReconstructedHitFactory.h index 0fcc5e4..aa69c9e 100644 --- a/source/tdis/tracking/ReconstructedHitFactory.h +++ b/source/tdis/tracking/ReconstructedHitFactory.h @@ -9,99 +9,89 @@ #include "podio_model/TrackerHit.h" #include "podio_model/TrackerHitCollection.h" -namespace tdis::tracking { - +namespace { + inline double get_resolution(const double pixel_size) { + constexpr const double sqrt_12 = 3.4641016151; + return pixel_size / sqrt_12; + } + inline double get_variance(const double pixel_size) { + const double res = get_resolution(pixel_size); + return res * res; + } +} // namespace +namespace tdis::tracking { struct ReconstructedHitFactory : public JOmniFactory { - - PodioInput m_mc_hits_in {this}; - PodioOutput m_tracker_hits_out {this}; + PodioInput m_mc_hits_in {this, {"DigitizedMtpcMcHit"}}; + PodioOutput m_tracker_hits_out {this, "TrackerHit"}; Service m_service_geometry{this}; - + Parameter m_cfg_use_true_pos{this, "acts:use_true_position", true,"Use true hits xyz instead of digitized one"}; void Configure() { m_service_geometry(); - } void ChangeRun(int32_t /*run_nr*/) { } void Execute(int32_t /*run_nr*/, uint64_t /*evt_nr*/) { + using namespace Acts::UnitLiterals; - auto result_hits = std::make_unique(); - - for (auto mc_hit : *m_mc_hits_in()) { - - auto [x,y] = getPadCenter(mc_hit.ring(), mc_hit.pad()); - //#double z = getRing - - //std::uint64_t cellID, edm4hep::Vector3f position, edm4eic::CovDiag3f positionError, float time, float timeError, float edep, float edepError - //cluster.addClusters(protocluster); - //cs->push_back(cluster); - } - - m_tracker_hits_out() = std::move(result_hits); - - /* - *auto rec_hits { std::make_unique() }; - - for (const auto& raw_hit : raw_hits) { - auto id = raw_hit.getCellID(); + auto rec_hits = std::make_unique(); + auto plane_positions = m_service_geometry->GetPlanePositions(); - // Get position and dimension - auto pos = m_converter->position(id); - auto dim = m_converter->cellDimensions(id); + for (auto mc_hit : *m_mc_hits_in()) { - // >oO trace - if(m_log->level() == spdlog::level::trace) { - m_log->trace("position x={:.2f} y={:.2f} z={:.2f} [mm]: ", pos.x()/ mm, pos.y()/ mm, pos.z()/ mm); - m_log->trace("dimension size: {}", dim.size()); - for (size_t j = 0; j < std::size(dim); ++j) { - m_log->trace(" - dimension {:<5} size: {:.2}", j, dim[j]); + const int plane = mc_hit.plane(); + const int ring = mc_hit.ring(); + const int pad = mc_hit.pad(); + const double z_to_gem = mc_hit.zToGem(); + + auto [pad_x,pad_y] = getPadCenter(ring, pad); + double plane_z = m_service_geometry().GetPlanePositions()[plane]; + + + // Planes are set back to back like this: + // | || || || | + // This means that zToPlane is positive for even planes and negative for odd + // i.e. + // z = plane_z + zToPlane // for even + // z = plane_z - zToPlane // for odd + double calc_z = plane_z + (plane % 2 ? -z_to_gem : z_to_gem); + + + // Position + edm4hep::Vector3f position; + if(m_cfg_use_true_pos() && !std::isnan(mc_hit.truePosition().x)) { + position.x = mc_hit.truePosition().x; + position.y = mc_hit.truePosition().y; + position.z = mc_hit.truePosition().z; + } + else { + position.x = pad_x; + position.y = pad_y; + position.z = calc_z; + } + + // Covariance + double max_dimension = std::max(getPadApproxWidth(ring), getPadHight()); + double xy_variance = get_variance(max_dimension); + edm4eic::CovDiag3f cov{xy_variance, xy_variance, 1_cm}; // TODO correct Z variance + + uint32_t cell_id = 1000000*mc_hit.plane() + 1000*mc_hit.ring() + mc_hit.pad(); + + rec_hits->create( + cell_id, + position, + cov, + static_cast(mc_hit.time()), + static_cast(1_ns), // TODO correct time resolution + static_cast(mc_hit.adc()), + 0.0F); + m_tracker_hits_out() = std::move(rec_hits); } } - - // Note about variance: - // The variance is used to obtain a diagonal covariance matrix. - // Note that the covariance matrix is written in DD4hep surface coordinates, - // *NOT* global position coordinates. This implies that: - // - XY segmentation: xx -> sigma_x, yy-> sigma_y, zz -> 0, tt -> 0 - // - XZ segmentation: xx -> sigma_x, yy-> sigma_z, zz -> 0, tt -> 0 - // - 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. - #if EDM4EIC_VERSION_MAJOR >= 7 - auto rec_hit = - #endif - rec_hits->create( - raw_hit.getCellID(), // Raw DD4hep cell ID - edm4hep::Vector3f{static_cast(pos.x() / mm), static_cast(pos.y() / mm), static_cast(pos.z() / mm)}, // mm - edm4eic::CovDiag3f{get_variance(dim[0] / mm), get_variance(dim[1] / mm), // variance (see note above) - std::size(dim) > 2 ? get_variance(dim[2] / mm) : 0.}, - static_cast((double)(raw_hit.getTimeStamp()) / 1000.0), // ns - m_cfg.timeResolution, // in ns - static_cast(raw_hit.getCharge() / 1.0e6), // Collected energy (GeV) - 0.0F); // Error on the energy -#if EDM4EIC_VERSION_MAJOR >= 7 - rec_hit.setRawHit(raw_hit); -#endif - * */ - } - - private: - static inline double get_resolution(const double pixel_size) { - constexpr const double sqrt_12 = 3.4641016151; - return pixel_size / sqrt_12; - } - - static inline double get_variance(const double pixel_size) { - const double res = get_resolution(pixel_size); - return res * res; - } - }; - } \ No newline at end of file From 3f809485f9b05d5e0c202cede1518a92df5064c8 Mon Sep 17 00:00:00 2001 From: Dmitry Romanov Date: Mon, 18 Nov 2024 16:06:34 -0500 Subject: [PATCH 6/8] Making tracking work --- source/tdis/CMakeLists.txt | 11 +- source/tdis/io/PodioWriteProcessor.hpp | 20 +-- source/tdis/tdis_main.cpp | 7 +- source/tdis/tracking/ActsGeometryService.cc | 15 +- source/tdis/tracking/ActsGeometryService.h | 6 +- ...copeDetector.cpp => BuildMtpcDetector.cpp} | 47 +++-- ...copeDetector.hpp => BuildMtpcDetector.hpp} | 32 ++-- source/tdis/tracking/MtpcDetectorElement.cpp | 27 +++ source/tdis/tracking/MtpcDetectorElement.hpp | 142 +++++++++++++++ .../tdis/tracking/ReconstructedHitFactory.h | 2 +- source/tdis/tracking/TelescopeDetector.cpp | 62 ------- source/tdis/tracking/TelescopeDetector.hpp | 64 ------- .../tracking/TelescopeDetectorElement.cpp | 42 ----- .../tracking/TelescopeDetectorElement.hpp | 162 ------------------ 14 files changed, 235 insertions(+), 404 deletions(-) rename source/tdis/tracking/{BuildTelescopeDetector.cpp => BuildMtpcDetector.cpp} (73%) rename source/tdis/tracking/{BuildTelescopeDetector.hpp => BuildMtpcDetector.hpp} (70%) create mode 100644 source/tdis/tracking/MtpcDetectorElement.cpp create mode 100644 source/tdis/tracking/MtpcDetectorElement.hpp delete mode 100644 source/tdis/tracking/TelescopeDetector.cpp delete mode 100644 source/tdis/tracking/TelescopeDetector.hpp delete mode 100644 source/tdis/tracking/TelescopeDetectorElement.cpp delete mode 100644 source/tdis/tracking/TelescopeDetectorElement.hpp diff --git a/source/tdis/CMakeLists.txt b/source/tdis/CMakeLists.txt index f6d3e66..859f509 100644 --- a/source/tdis/CMakeLists.txt +++ b/source/tdis/CMakeLists.txt @@ -43,13 +43,10 @@ add_executable(tdis tracking/ActsGeometryService.cc tracking/ActsGeometryService.h tracking/ReconstructedHitFactory.h - - tracking/BuildTelescopeDetector.cpp - tracking/BuildTelescopeDetector.hpp - tracking/TelescopeDetector.cpp - tracking/TelescopeDetector.hpp - tracking/TelescopeDetectorElement.cpp - tracking/TelescopeDetectorElement.hpp + tracking/BuildMtpcDetector.cpp + tracking/BuildMtpcDetector.hpp + tracking/MtpcDetectorElement.cpp + tracking/MtpcDetectorElement.hpp # tracking/CKFTrackingFunction.cc # tracking/DD4hepBField.h # /tracking/DD4hepBField.cc diff --git a/source/tdis/io/PodioWriteProcessor.hpp b/source/tdis/io/PodioWriteProcessor.hpp index 690a4be..1278d51 100644 --- a/source/tdis/io/PodioWriteProcessor.hpp +++ b/source/tdis/io/PodioWriteProcessor.hpp @@ -36,19 +36,9 @@ * repeated exceptions and potential crashes. */ +#include #include #include -#include -#include -#include -#include -#include -#include -#include -#include - - -#include #include #include #include @@ -57,11 +47,18 @@ #include #include #include +#include #include #include +#include +#include +#include +#include #include +#include +#include "podio_model/TrackerHit.h" #include "services/LogService.hpp" namespace tdis::io { @@ -142,6 +139,7 @@ inline void PodioWriteProcessor::Init() { inline void PodioWriteProcessor::Process(const std::shared_ptr& event) { std::lock_guard lock(m_mutex); + auto hits = event->GetCollection("TrackerHit"); m_log->info("PodioWriteProcessor::Process() All event collections:"); auto event_collections = event->GetAllCollectionNames(); diff --git a/source/tdis/tdis_main.cpp b/source/tdis/tdis_main.cpp index 601e599..d40fe55 100644 --- a/source/tdis/tdis_main.cpp +++ b/source/tdis/tdis_main.cpp @@ -11,6 +11,7 @@ #include #include +#include #include @@ -114,9 +115,13 @@ int main(int argc, char* argv[]) { app.ProvideService(std::make_shared(&app)); app.ProvideService(std::make_shared()); + auto reco_hit_generator = new JOmniFactoryGeneratorT(); + reco_hit_generator->AddWiring("TrackerHitGenerator", {"DigitizedMtpcMcHit"}, {"TrackerHit"}); + app.Add(reco_hit_generator); + app.Add(new JEventSourceGeneratorT); app.Add(new tdis::io::PodioWriteProcessor(&app)); - app.Add(new JFactoryGeneratorT("TrackerHit")); + // app.Add(new JEventProcessorPodio); // app.Add(new JFactoryGeneratorT()); // app.Add(new JFactoryGeneratorT()); diff --git a/source/tdis/tracking/ActsGeometryService.cc b/source/tdis/tracking/ActsGeometryService.cc index 8000952..f6fb25d 100644 --- a/source/tdis/tracking/ActsGeometryService.cc +++ b/source/tdis/tracking/ActsGeometryService.cc @@ -23,8 +23,7 @@ #include #include -#include "BuildTelescopeDetector.hpp" -#include "TelescopeDetector.hpp" +#include "BuildMtpcDetector.hpp" #include "services/LogService.hpp" // Formatter for Eigen matrices @@ -369,19 +368,19 @@ void tdis::tracking::ActsGeometryService::Init() { /// Return the telescope detector - gGeometry = - ActsExamples::Telescope::buildDetector( + gGeometry = tdis::tracking::buildDetector( nominalContext, // gctx is the detector element dependent geometry context - detectorStore, // detectorStore is the store for the detector element - m_plane_positions, // positions are the positions of different layers in the longitudinal direction + detectorStore, // detectorStore is the store for the detector element + m_plane_positions, // positions are the positions of different layers in the longitudinal direction stereos, // stereoAngles are the stereo angles of different layers, which are rotation angles around the longitudinal (normal) direction offsets, // is the offset (u, v) of the layers in the transverse plane - bounds, // bounds is the surface bound values, i.e. halfX and halfY if plane surface, and minR and maxR if disc surface + bounds, // bounds is the surface bound values, minR and maxR thickness, // thickness is the material thickness of each layer - ActsExamples::Telescope::TelescopeSurfaceType::Disc, // surfaceType is the detector surface type Acts::BinningValue::binZ); + + // Visualize ACTS geometry const Acts::TrackingVolume& tgVolume = *(gGeometry->highestTrackingVolume()); diff --git a/source/tdis/tracking/ActsGeometryService.h b/source/tdis/tracking/ActsGeometryService.h index fe99e27..c6b9367 100644 --- a/source/tdis/tracking/ActsGeometryService.h +++ b/source/tdis/tracking/ActsGeometryService.h @@ -14,7 +14,7 @@ #include #include -#include "TelescopeDetectorElement.hpp" +#include "MtpcDetectorElement.hpp" #include "services/LogService.hpp" namespace tdis::tracking { @@ -60,9 +60,9 @@ namespace tdis::tracking { Acts::GeometryContext m_geometry_context = Acts::GeometryContext(); - ActsExamples::Telescope::TelescopeDetectorElement::ContextType nominalContext; + tdis::tracking::MtpcDetectorElement::ContextType nominalContext; - std::vector> + std::vector> detectorStore; std::shared_ptr gGeometry; diff --git a/source/tdis/tracking/BuildTelescopeDetector.cpp b/source/tdis/tracking/BuildMtpcDetector.cpp similarity index 73% rename from source/tdis/tracking/BuildTelescopeDetector.cpp rename to source/tdis/tracking/BuildMtpcDetector.cpp index 3f7835c..4e7af2e 100644 --- a/source/tdis/tracking/BuildTelescopeDetector.cpp +++ b/source/tdis/tracking/BuildMtpcDetector.cpp @@ -6,7 +6,9 @@ // License, v. 2.0. If a copy of the MPL was not distributed with this // file, You can obtain one at http://mozilla.org/MPL/2.0/. -#include "BuildTelescopeDetector.hpp" +#include +#include +#include #include "Acts/Definitions/Algebra.hpp" #include "Acts/Definitions/Units.hpp" @@ -28,17 +30,22 @@ #include "Acts/Surfaces/Surface.hpp" #include "Acts/Surfaces/SurfaceArray.hpp" #include "Acts/Utilities/Logger.hpp" - -#include -#include -#include +#include "BuildMtpcDetector.hpp" std::unique_ptr -ActsExamples::Telescope::buildDetector(const typename ActsExamples::Telescope::TelescopeDetectorElement::ContextType &gctx, std::vector > &detectorStore, const std::vector &positions, const std::vector &stereoAngles, const std::array &offsets, const std::array &bounds, double thickness, ActsExamples::Telescope::TelescopeSurfaceType surfaceType, Acts::BinningValue binValue) { +tdis::tracking::buildDetector( + const typename tdis::tracking::MtpcDetectorElement::ContextType &gctx, + std::vector > &detectorStore, + const std::vector &positions, + const std::vector &stereoAngles, + const std::array &offsets, + const std::array &bounds, + double thickness, + Acts::BinningValue binValue) +{ using namespace Acts::UnitLiterals; - // The rectangle bounds for plane surface - const auto pBounds = std::make_shared(bounds[0], bounds[1]); + // The radial bounds for disc surface const auto rBounds = std::make_shared(bounds[0], bounds[1]); @@ -76,24 +83,20 @@ ActsExamples::Telescope::buildDetector(const typename ActsExamples::Telescope::T trafo *= Acts::AngleAxis3(stereo, Acts::Vector3::UnitZ()); // Create the detector element - std::shared_ptr detElement = nullptr; - if (surfaceType == TelescopeSurfaceType::Plane) { - detElement = std::make_shared(std::make_shared(trafo), pBounds, 1._um, surfaceMaterial); - } else { - detElement = std::make_shared(std::make_shared(trafo), rBounds, 1._um, surfaceMaterial); - } + std::shared_ptr detElement = nullptr; + + detElement = std::make_shared(std::make_shared(trafo), rBounds, 1._um, i, surfaceMaterial); + // Get the surface auto surface = detElement->surface().getSharedPtr(); + // Add the detector element to the detector store detectorStore.push_back(std::move(detElement)); // Construct the surface array (one surface contained) std::unique_ptr surArray(new Acts::SurfaceArray(surface)); // Construct the layer - if (surfaceType == TelescopeSurfaceType::Plane) { - layers[i] = Acts::PlaneLayer::create(trafo, pBounds, std::move(surArray), thickness); - } else { - layers[i] = Acts::DiscLayer::create(trafo, rBounds, std::move(surArray), thickness); - } + layers[i] = Acts::DiscLayer::create(trafo, rBounds, std::move(surArray), thickness); + // Associate the layer to the surface auto mutableSurface = const_cast(surface.get()); mutableSurface->associateLayer(*layers[i]); @@ -107,11 +110,7 @@ ActsExamples::Telescope::buildDetector(const typename ActsExamples::Telescope::T // or cylinder with discs auto length = positions.back() - positions.front(); std::shared_ptr boundsVol = nullptr; - if (surfaceType == TelescopeSurfaceType::Plane) { - boundsVol = std::make_shared(bounds[0] + 5._mm, bounds[1] + 5._mm, length + 10._mm); - } else { - boundsVol = std::make_shared(std::max(bounds[0] - 5.0_mm, 0.), bounds[1] + 5._mm, length + 10._mm); - } + boundsVol = std::make_shared(std::max(bounds[0] - 5.0_mm, 0.), bounds[1] + 5._mm, length + 10._mm); Acts::LayerArrayCreator::Config lacConfig; Acts::LayerArrayCreator layArrCreator(lacConfig, Acts::getDefaultLogger("LayerArrayCreator", Acts::Logging::INFO)); diff --git a/source/tdis/tracking/BuildTelescopeDetector.hpp b/source/tdis/tracking/BuildMtpcDetector.hpp similarity index 70% rename from source/tdis/tracking/BuildTelescopeDetector.hpp rename to source/tdis/tracking/BuildMtpcDetector.hpp index f64c5c3..5f8cb7e 100644 --- a/source/tdis/tracking/BuildTelescopeDetector.hpp +++ b/source/tdis/tracking/BuildMtpcDetector.hpp @@ -8,26 +8,20 @@ #pragma once -#include "Acts/Geometry/TrackingGeometry.hpp" -#include "Acts/Utilities/BinUtility.hpp" -#include "Acts/Utilities/BinningType.hpp" -#include "TelescopeDetectorElement.hpp" - #include #include #include +#include "Acts/Geometry/TrackingGeometry.hpp" +#include "Acts/Utilities/BinUtility.hpp" +#include "Acts/Utilities/BinningType.hpp" +#include "MtpcDetectorElement.hpp" + namespace Acts { class TrackingGeometry; } // namespace Acts -namespace ActsExamples::Telescope { - -/// The telescope detector surface type -enum class TelescopeSurfaceType { - Plane = 0, - Disc = 1, -}; +namespace tdis::tracking { /// Global method to build the telescope tracking geometry /// @@ -39,19 +33,19 @@ enum class TelescopeSurfaceType { /// rotation angles around the longitudinal (normal) /// direction /// @param offsets is the offset (u, v) of the layers in the transverse plane -/// @param bounds is the surface bound values, i.e. halfX and halfY if plane -/// surface, and minR and maxR if disc surface +/// @param bounds is the surface bound values, i.e. minR and maxR /// @param thickness is the material thickness of each layer -/// @param surfaceType is the detector surface type + /// @param binValue indicates which axis the detector surface normals are /// parallel to std::unique_ptr buildDetector( - const typename TelescopeDetectorElement::ContextType& gctx, - std::vector>& detectorStore, + const typename MtpcDetectorElement::ContextType& gctx, + std::vector>& detectorStore, const std::vector& positions, const std::vector& stereoAngles, - const std::array& offsets, const std::array& bounds, - double thickness, TelescopeSurfaceType surfaceType, + const std::array& offsets, + const std::array& bounds, + double thickness, Acts::BinningValue binValue = Acts::BinningValue::binZ); } // namespace ActsExamples::Telescope diff --git a/source/tdis/tracking/MtpcDetectorElement.cpp b/source/tdis/tracking/MtpcDetectorElement.cpp new file mode 100644 index 0000000..e256de3 --- /dev/null +++ b/source/tdis/tracking/MtpcDetectorElement.cpp @@ -0,0 +1,27 @@ +// This file is part of the Acts project. +// +// Copyright (C) 2020-2021 CERN for the benefit of the Acts project +// +// This Source Code Form is subject to the terms of the Mozilla Public +// License, v. 2.0. If a copy of the MPL was not distributed with this +// file, You can obtain one at http://mozilla.org/MPL/2.0/. + +#include "Acts/Surfaces/DiscSurface.hpp" +#include "Acts/Surfaces/PlaneSurface.hpp" +#include "MtpcDetectorElement.hpp" + + +tdis::tracking::MtpcDetectorElement::MtpcDetectorElement( + std::shared_ptr transform, + std::shared_ptr dBounds, + double thickness, + int gem_plane_id, + std::shared_ptr material) + : Acts::DetectorElementBase(), + m_elementTransform(std::move(transform)), + m_elementSurface(Acts::Surface::makeShared(dBounds, *this)), + m_elementThickness(thickness), + m_gem_plane_id(gem_plane_id), + m_elementDiscBounds(std::move(dBounds)) { + m_elementSurface->assignSurfaceMaterial(std::move(material)); +} diff --git a/source/tdis/tracking/MtpcDetectorElement.hpp b/source/tdis/tracking/MtpcDetectorElement.hpp new file mode 100644 index 0000000..eb955ef --- /dev/null +++ b/source/tdis/tracking/MtpcDetectorElement.hpp @@ -0,0 +1,142 @@ +// This file is part of the Acts project. +// +// Copyright (C) 2020-2021 CERN for the benefit of the Acts project +// +// This Source Code Form is subject to the terms of the Mozilla Public +// License, v. 2.0. If a copy of the MPL was not distributed with this +// file, You can obtain one at http://mozilla.org/MPL/2.0/. + +#pragma once + +#include +#include +#include + +#include "Acts/Definitions/Algebra.hpp" +#include "Acts/Geometry/DetectorElementBase.hpp" +#include "Acts/Geometry/GeometryContext.hpp" +#include "Acts/Surfaces/Surface.hpp" + +namespace Acts { + class Surface; + class PlanarBounds; + class DiscBounds; + class ISurfaceMaterial; +} // namespace Acts + +namespace tdis::tracking { + + /// @class MtpcDetectorElement + /// + /// This is a lightweight type of detector element, + /// it simply implements the base class. + class MtpcDetectorElement : public Acts::DetectorElementBase { + public: + /// @class ContextType + /// convention: nested to the Detector element + struct ContextType { + /// The current interval of validity + unsigned int iov = 0; + }; + + /// Constructor for single sided detector element + /// - bound to a Disc Surface + /// + /// @param transform is the transform that element the layer in 3D frame + /// @param dBounds is the planar bounds for the disc like detector element + /// @param thickness is the module thickness + /// @param gem_plane_id the id of the GEM plane + /// @param material is the (optional) Surface material associated to it + MtpcDetectorElement(std::shared_ptr transform, + std::shared_ptr dBounds, double thickness, + int gem_plane_id, + std::shared_ptr material = nullptr + ); + + /// Destructor + ~MtpcDetectorElement() override = default; + + /// Return surface associated with this detector element + const Acts::Surface& surface() const final; + + /// Non-const access to the surface associated with this detector element + Acts::Surface& surface() final; + + /// The maximal thickness of the detector element wrt normal axis + double thickness() const final; + + /// Return local to global transform associated with this identifier + /// + /// @param gctx The current geometry context object, e.g. alignment + /// + /// @note this is called from the surface().transform() in the PROXY mode + const Acts::Transform3& transform(const Acts::GeometryContext& gctx) const final; + + /// Return the nominal local to global transform + /// + /// @note the geometry context will hereby be ignored + const Acts::Transform3& nominalTransform(const Acts::GeometryContext& gctx) const; + + /// Return local to global transform associated with this identifier + /// + /// @param alignedTransform is a new transform + /// @param iov is the batch for which it is meant + void addAlignedTransform(std::unique_ptr alignedTransform, unsigned int iov); + + /// Return the set of alignment transforms in flight + const std::vector>& alignedTransforms() const; + + private: + /// the transform for positioning in 3D space + std::shared_ptr m_elementTransform = nullptr; + // the aligned transforms + std::vector> m_alignedTransforms = {}; + /// the surface represented by it + std::shared_ptr m_elementSurface = nullptr; + /// the element thickness + double m_elementThickness = 0.; + /// the disc bounds + std::shared_ptr m_elementDiscBounds = nullptr; + /// GEM plane id + int m_gem_plane_id; + }; + + inline const Acts::Surface& MtpcDetectorElement::surface() const { return *m_elementSurface; } + + inline Acts::Surface& MtpcDetectorElement::surface() { return *m_elementSurface; } + + inline double MtpcDetectorElement::thickness() const { return m_elementThickness; } + + inline const Acts::Transform3& MtpcDetectorElement::transform(const Acts::GeometryContext& gctx) const + { + // Check if a different transform than the nominal exists + if (!m_alignedTransforms.empty()) { + // cast into the right context object + auto alignContext = gctx.get(); + return (*m_alignedTransforms[alignContext.iov].get()); + } + // Return the standard transform if not found + return nominalTransform(gctx); + } + + inline const Acts::Transform3& MtpcDetectorElement::nominalTransform(const Acts::GeometryContext& /*gctx*/) const + { + return *m_elementTransform; + } + + inline void MtpcDetectorElement::addAlignedTransform(std::unique_ptr alignedTransform, unsigned int iov) + { + // most standard case, it's just a new one: + auto cios = m_alignedTransforms.size(); + for (unsigned int ic = cios; ic <= iov; ++ic) { + m_alignedTransforms.push_back(nullptr); + } + m_alignedTransforms[iov] = std::move(alignedTransform); + } + + inline const std::vector>& + MtpcDetectorElement::alignedTransforms() const { + return m_alignedTransforms; + } + +} // namespace tdis::tracking diff --git a/source/tdis/tracking/ReconstructedHitFactory.h b/source/tdis/tracking/ReconstructedHitFactory.h index aa69c9e..15cf2e5 100644 --- a/source/tdis/tracking/ReconstructedHitFactory.h +++ b/source/tdis/tracking/ReconstructedHitFactory.h @@ -90,8 +90,8 @@ namespace tdis::tracking { static_cast(1_ns), // TODO correct time resolution static_cast(mc_hit.adc()), 0.0F); - m_tracker_hits_out() = std::move(rec_hits); } + m_tracker_hits_out() = std::move(rec_hits); } }; } \ No newline at end of file diff --git a/source/tdis/tracking/TelescopeDetector.cpp b/source/tdis/tracking/TelescopeDetector.cpp deleted file mode 100644 index ba4f7bf..0000000 --- a/source/tdis/tracking/TelescopeDetector.cpp +++ /dev/null @@ -1,62 +0,0 @@ -// This file is part of the Acts project. -// -// Copyright (C) 2020-2021 CERN for the benefit of the Acts project -// -// This Source Code Form is subject to the terms of the Mozilla Public -// License, v. 2.0. If a copy of the MPL was not distributed with this -// file, You can obtain one at http://mozilla.org/MPL/2.0/. - -#include "TelescopeDetector.hpp" - -#include "Acts/Geometry/TrackingGeometry.hpp" -#include "Acts/Utilities/BinningType.hpp" -#include "BuildTelescopeDetector.hpp" -#include "TelescopeDetectorElement.hpp" - -#include -#include - -auto ActsExamples::Telescope::TelescopeDetector::finalize( - const Config& cfg, const std::shared_ptr& - /*mdecorator*/) -> std::pair { - DetectorElement::ContextType nominalContext; - - if (cfg.surfaceType > 1) { - throw std::invalid_argument( - "The surface type could either be 0 for plane surface or 1 for disc " - "surface."); - } - if (cfg.binValue > 2) { - throw std::invalid_argument("The axis value could only be 0, 1, or 2."); - } - // Check if the bounds values are valid - if (cfg.surfaceType == 1 && cfg.bounds[0] >= cfg.bounds[1]) { - throw std::invalid_argument( - "The minR should be smaller than the maxR for disc surface bounds."); - } - - if (cfg.positions.size() != cfg.stereos.size()) { - throw std::invalid_argument( - "The number of provided positions must match the number of " - "provided stereo angles."); - } - - config = cfg; - - // Sort the provided distances - std::vector positions = cfg.positions; - std::vector stereos = cfg.stereos; - std::sort(positions.begin(), positions.end()); - - /// Return the telescope detector - TrackingGeometryPtr gGeometry = ActsExamples::Telescope::buildDetector( - nominalContext, detectorStore, positions, stereos, cfg.offsets, - cfg.bounds, cfg.thickness, - static_cast( - cfg.surfaceType), - static_cast(cfg.binValue)); - ContextDecorators gContextDecorators = {}; - // return the pair of geometry and empty decorators - return std::make_pair( - std::move(gGeometry), std::move(gContextDecorators)); -} diff --git a/source/tdis/tracking/TelescopeDetector.hpp b/source/tdis/tracking/TelescopeDetector.hpp deleted file mode 100644 index 19bb2b7..0000000 --- a/source/tdis/tracking/TelescopeDetector.hpp +++ /dev/null @@ -1,64 +0,0 @@ -// This file is part of the Acts project. -// -// Copyright (C) 2020 CERN for the benefit of the Acts project -// -// This Source Code Form is subject to the terms of the Mozilla Public -// License, v. 2.0. If a copy of the MPL was not distributed with this -// file, You can obtain one at http://mozilla.org/MPL/2.0/. - -#pragma once - -#include "Acts/Definitions/Units.hpp" -#include "Acts/Utilities/Logger.hpp" -#include "ActsExamples/Utilities/Options.hpp" - -#include -#include -#include -#include - -using namespace Acts::UnitLiterals; - -namespace Acts { -class TrackingGeometry; -class IMaterialDecorator; -} // namespace Acts - -namespace ActsExamples { -class IContextDecorator; -} // namespace ActsExamples - -namespace ActsExamples::Telescope { - -class TelescopeDetectorElement; -class TelescopeG4DetectorConstruction; - -struct TelescopeDetector { - using DetectorElement = ActsExamples::Telescope::TelescopeDetectorElement; - using DetectorElementPtr = std::shared_ptr; - using DetectorStore = std::vector; - - using ContextDecorators = - std::vector>; - using TrackingGeometryPtr = std::shared_ptr; - - struct Config { - std::vector positions{{0, 30, 60, 120, 150, 180}}; - std::vector stereos{{0, 0, 0, 0, 0, 0}}; - std::array offsets{{0, 0}}; - std::array bounds{{25, 100}}; - double thickness{80_um}; - int surfaceType{0}; - int binValue{2}; - }; - - Config config; - /// The store of the detector elements (lifetime: job) - DetectorStore detectorStore; - - std::pair finalize( - const Config& cfg, - const std::shared_ptr& mdecorator); -}; - -} // namespace ActsExamples::Telescope diff --git a/source/tdis/tracking/TelescopeDetectorElement.cpp b/source/tdis/tracking/TelescopeDetectorElement.cpp deleted file mode 100644 index c761d43..0000000 --- a/source/tdis/tracking/TelescopeDetectorElement.cpp +++ /dev/null @@ -1,42 +0,0 @@ -// This file is part of the Acts project. -// -// Copyright (C) 2020-2021 CERN for the benefit of the Acts project -// -// This Source Code Form is subject to the terms of the Mozilla Public -// License, v. 2.0. If a copy of the MPL was not distributed with this -// file, You can obtain one at http://mozilla.org/MPL/2.0/. - -#include "TelescopeDetectorElement.hpp" - -#include "Acts/Surfaces/DiscSurface.hpp" -#include "Acts/Surfaces/PlaneSurface.hpp" - -ActsExamples::Telescope::TelescopeDetectorElement::TelescopeDetectorElement( - std::shared_ptr transform, - std::shared_ptr pBounds, double thickness, - std::shared_ptr material) - : Acts::DetectorElementBase(), - m_elementTransform(std::move(transform)), - m_elementSurface( - Acts::Surface::makeShared(pBounds, *this)), - m_elementThickness(thickness), - m_elementPlanarBounds(std::move(pBounds)), - m_elementDiscBounds(nullptr) { - auto mutableSurface = - std::const_pointer_cast(m_elementSurface); - mutableSurface->assignSurfaceMaterial(std::move(material)); -} - -ActsExamples::Telescope::TelescopeDetectorElement::TelescopeDetectorElement( - std::shared_ptr transform, - std::shared_ptr dBounds, double thickness, - std::shared_ptr material) - : Acts::DetectorElementBase(), - m_elementTransform(std::move(transform)), - m_elementSurface( - Acts::Surface::makeShared(dBounds, *this)), - m_elementThickness(thickness), - m_elementPlanarBounds(nullptr), - m_elementDiscBounds(std::move(dBounds)) { - m_elementSurface->assignSurfaceMaterial(std::move(material)); -} diff --git a/source/tdis/tracking/TelescopeDetectorElement.hpp b/source/tdis/tracking/TelescopeDetectorElement.hpp deleted file mode 100644 index e56bbf9..0000000 --- a/source/tdis/tracking/TelescopeDetectorElement.hpp +++ /dev/null @@ -1,162 +0,0 @@ -// This file is part of the Acts project. -// -// Copyright (C) 2020-2021 CERN for the benefit of the Acts project -// -// This Source Code Form is subject to the terms of the Mozilla Public -// License, v. 2.0. If a copy of the MPL was not distributed with this -// file, You can obtain one at http://mozilla.org/MPL/2.0/. - -#pragma once - -#include "Acts/Definitions/Algebra.hpp" -#include "Acts/Geometry/DetectorElementBase.hpp" -#include "Acts/Geometry/GeometryContext.hpp" -#include "Acts/Surfaces/Surface.hpp" - -#include -#include -#include - -namespace Acts { -class Surface; -class PlanarBounds; -class DiscBounds; -class ISurfaceMaterial; -} // namespace Acts - -namespace ActsExamples::Telescope { - -/// @class TelescopeDetectorElement -/// -/// This is a lightweight type of detector element, -/// it simply implements the base class. -class TelescopeDetectorElement : public Acts::DetectorElementBase { - public: - /// @class ContextType - /// convention: nested to the Detector element - struct ContextType { - /// The current interval of validity - unsigned int iov = 0; - }; - - /// Constructor for single sided detector element - /// - bound to a Plane Surface - /// - /// @param transform is the transform that element the layer in 3D frame - /// @param pBounds is the planar bounds for the planar detector element - /// @param thickness is the module thickness - /// @param material is the (optional) Surface material associated to it - TelescopeDetectorElement( - std::shared_ptr transform, - std::shared_ptr pBounds, double thickness, - std::shared_ptr material = nullptr); - - /// Constructor for single sided detector element - /// - bound to a Disc Surface - /// - /// @param transform is the transform that element the layer in 3D frame - /// @param dBounds is the planar bounds for the disc like detector element - /// @param thickness is the module thickness - /// @param material is the (optional) Surface material associated to it - TelescopeDetectorElement( - std::shared_ptr transform, - std::shared_ptr dBounds, double thickness, - std::shared_ptr material = nullptr); - - /// Destructor - ~TelescopeDetectorElement() override = default; - - /// Return surface associated with this detector element - const Acts::Surface& surface() const final; - - /// Non-const access to the surface associated with this detector element - Acts::Surface& surface() final; - - /// The maximal thickness of the detector element wrt normal axis - double thickness() const final; - - /// Return local to global transform associated with this identifier - /// - /// @param gctx The current geometry context object, e.g. alignment - /// - /// @note this is called from the surface().transform() in the PROXY mode - const Acts::Transform3& transform( - const Acts::GeometryContext& gctx) const final; - - /// Return the nominal local to global transform - /// - /// @note the geometry context will hereby be ignored - const Acts::Transform3& nominalTransform( - const Acts::GeometryContext& gctx) const; - - /// Return local to global transform associated with this identifier - /// - /// @param alignedTransform is a new transform - /// @param iov is the batch for which it is meant - void addAlignedTransform(std::unique_ptr alignedTransform, - unsigned int iov); - - /// Return the set of alignment transforms in flight - const std::vector>& alignedTransforms() - const; - - private: - /// the transform for positioning in 3D space - std::shared_ptr m_elementTransform = nullptr; - // the aligned transforms - std::vector> m_alignedTransforms = {}; - /// the surface represented by it - std::shared_ptr m_elementSurface = nullptr; - /// the element thickness - double m_elementThickness = 0.; - /// the planar bounds - std::shared_ptr m_elementPlanarBounds = nullptr; - /// the disc bounds - std::shared_ptr m_elementDiscBounds = nullptr; -}; - -inline const Acts::Surface& TelescopeDetectorElement::surface() const { - return *m_elementSurface; -} - -inline Acts::Surface& TelescopeDetectorElement::surface() { - return *m_elementSurface; -} - -inline double TelescopeDetectorElement::thickness() const { - return m_elementThickness; -} - -inline const Acts::Transform3& TelescopeDetectorElement::transform( - const Acts::GeometryContext& gctx) const { - // Check if a different transform than the nominal exists - if (!m_alignedTransforms.empty()) { - // cast into the right context object - auto alignContext = gctx.get(); - return (*m_alignedTransforms[alignContext.iov].get()); - } - // Return the standard transform if not found - return nominalTransform(gctx); -} - -inline const Acts::Transform3& TelescopeDetectorElement::nominalTransform( - const Acts::GeometryContext& /*gctx*/) const { - return *m_elementTransform; -} - -inline void TelescopeDetectorElement::addAlignedTransform( - std::unique_ptr alignedTransform, unsigned int iov) { - // most standard case, it's just a new one: - auto cios = m_alignedTransforms.size(); - for (unsigned int ic = cios; ic <= iov; ++ic) { - m_alignedTransforms.push_back(nullptr); - } - m_alignedTransforms[iov] = std::move(alignedTransform); -} - -inline const std::vector>& -TelescopeDetectorElement::alignedTransforms() const { - return m_alignedTransforms; -} - -} // namespace ActsExamples::Telescope From 2a501db4a9514345dfbb6eab32c34987df4c4b0f Mon Sep 17 00:00:00 2001 From: Dmitry Romanov Date: Tue, 19 Nov 2024 09:37:19 -0500 Subject: [PATCH 7/8] 2D Measurements --- source/tdis/CMakeLists.txt | 7 +- source/tdis/tdis_main.cpp | 10 +- source/tdis/tracking/ActsGeometryService.cc | 2 +- source/tdis/tracking/ActsGeometryService.h | 8 +- source/tdis/tracking/CkfFitFactory.h | 36 ------ source/tdis/tracking/KalmanFittingFactory.h | 24 ++++ source/tdis/tracking/Measurement2DFactory.h | 41 +++++++ .../tdis/tracking/ReconstructedHitFactory.h | 111 ++++++++++++++---- .../tracking/TruthTrackParameterFactory.h | 97 +++++++++++++++ 9 files changed, 268 insertions(+), 68 deletions(-) delete mode 100644 source/tdis/tracking/CkfFitFactory.h create mode 100644 source/tdis/tracking/KalmanFittingFactory.h create mode 100644 source/tdis/tracking/Measurement2DFactory.h create mode 100644 source/tdis/tracking/TruthTrackParameterFactory.h diff --git a/source/tdis/CMakeLists.txt b/source/tdis/CMakeLists.txt index 859f509..2c17dec 100644 --- a/source/tdis/CMakeLists.txt +++ b/source/tdis/CMakeLists.txt @@ -47,9 +47,10 @@ add_executable(tdis tracking/BuildMtpcDetector.hpp tracking/MtpcDetectorElement.cpp tracking/MtpcDetectorElement.hpp - # tracking/CKFTrackingFunction.cc - # tracking/DD4hepBField.h - # /tracking/DD4hepBField.cc + # tracking/Measurement2DFactory.h + tracking/TruthTrackParameterFactory.h + tracking/KalmanFittingFactory.h + ) # ---------- FIND REQUIRED PACKAGES ------------- diff --git a/source/tdis/tdis_main.cpp b/source/tdis/tdis_main.cpp index d40fe55..0937994 100644 --- a/source/tdis/tdis_main.cpp +++ b/source/tdis/tdis_main.cpp @@ -21,6 +21,7 @@ #include "services/LogService.hpp" #include "tracking/ActsGeometryService.h" #include "tracking/ReconstructedHitFactory.h" +// #include "tracking/Measurement2DFactory.h" struct ProgramArguments { std::map params; @@ -116,9 +117,16 @@ int main(int argc, char* argv[]) { app.ProvideService(std::make_shared()); auto reco_hit_generator = new JOmniFactoryGeneratorT(); - reco_hit_generator->AddWiring("TrackerHitGenerator", {"DigitizedMtpcMcHit"}, {"TrackerHit"}); + reco_hit_generator->AddWiring( + "TrackerHitGenerator", + {"DigitizedMtpcMcHit"}, + {"TrackerHit", "Measurement2D"}); app.Add(reco_hit_generator); + // auto measurement_2d_generator = new JOmniFactoryGeneratorT(); + // measurement_2d_generator->AddWiring("TrackerHitGenerator", {"TrackerHit"}, {"Measurement2D"}); + // app.Add(measurement_2d_generator); + app.Add(new JEventSourceGeneratorT); app.Add(new tdis::io::PodioWriteProcessor(&app)); diff --git a/source/tdis/tracking/ActsGeometryService.cc b/source/tdis/tracking/ActsGeometryService.cc index f6fb25d..2e901f7 100644 --- a/source/tdis/tracking/ActsGeometryService.cc +++ b/source/tdis/tracking/ActsGeometryService.cc @@ -370,7 +370,7 @@ void tdis::tracking::ActsGeometryService::Init() { /// Return the telescope detector gGeometry = tdis::tracking::buildDetector( nominalContext, // gctx is the detector element dependent geometry context - detectorStore, // detectorStore is the store for the detector element + m_detector_elements, // detectorStore is the store for the detector element m_plane_positions, // positions are the positions of different layers in the longitudinal direction stereos, // stereoAngles are the stereo angles of different layers, which are rotation angles around the longitudinal (normal) direction offsets, // is the offset (u, v) of the layers in the transverse plane diff --git a/source/tdis/tracking/ActsGeometryService.h b/source/tdis/tracking/ActsGeometryService.h index c6b9367..24ea401 100644 --- a/source/tdis/tracking/ActsGeometryService.h +++ b/source/tdis/tracking/ActsGeometryService.h @@ -39,6 +39,11 @@ namespace tdis::tracking { const Acts::TrackingVolume& tVolume, const Acts::GeometryContext& gctx); + /// Basically returns GEM plane information + std::shared_ptr GetDetectorElement(size_t index) const { + return m_detector_elements[index]; + } + private: Parameter m_tgeo_file{this, "acts:geometry", "g4sbs_mtpc.root","TGeo filename with geometry for ACTS"}; Parameter m_material_map_file{this, "acts:material_map", "", "JSON/CBOR material map file path"}; @@ -62,8 +67,7 @@ namespace tdis::tracking { tdis::tracking::MtpcDetectorElement::ContextType nominalContext; - std::vector> - detectorStore; + std::vector> m_detector_elements; std::shared_ptr gGeometry; diff --git a/source/tdis/tracking/CkfFitFactory.h b/source/tdis/tracking/CkfFitFactory.h deleted file mode 100644 index f6c99a6..0000000 --- a/source/tdis/tracking/CkfFitFactory.h +++ /dev/null @@ -1,36 +0,0 @@ -#pragma once - -#include -#include - -namespace tdis::tracking { - - - - struct MyClusterFactory : public JOmniFactory { - - PodioInput m_protoclusters_in {this}; - PodioOutput m_clusters_out {this}; - - - void Configure() { - } - - void ChangeRun(int32_t /*run_nr*/) { - } - - void Execute(int32_t /*run_nr*/, uint64_t /*evt_nr*/) { - - auto cs = std::make_unique(); - - for (auto protocluster : *m_protoclusters_in()) { - auto cluster = MutableExampleCluster(protocluster.energy() + 1000); - cluster.addClusters(protocluster); - cs->push_back(cluster); - } - - m_clusters_out() = std::move(cs); - } - }; - -} \ No newline at end of file diff --git a/source/tdis/tracking/KalmanFittingFactory.h b/source/tdis/tracking/KalmanFittingFactory.h new file mode 100644 index 0000000..1069961 --- /dev/null +++ b/source/tdis/tracking/KalmanFittingFactory.h @@ -0,0 +1,24 @@ +#pragma once + +#include +#include + +namespace tdis::tracking { + + + + struct KalmanFittingFactory : public JOmniFactory { + + void Configure() { + } + + void ChangeRun(int32_t /*run_nr*/) { + } + + void Execute(int32_t /*run_nr*/, uint64_t /*evt_nr*/) { + + + } + }; + +} \ No newline at end of file diff --git a/source/tdis/tracking/Measurement2DFactory.h b/source/tdis/tracking/Measurement2DFactory.h new file mode 100644 index 0000000..3a57b4e --- /dev/null +++ b/source/tdis/tracking/Measurement2DFactory.h @@ -0,0 +1,41 @@ +#pragma once + +#include +#include + +#include "PadGeometryHelper.hpp" + +#include "podio_model/MutableTrackerHit.h" +#include "podio_model/TrackerHit.h" +#include "podio_model/TrackerHitCollection.h" +#include "podio_model/Measurement2D.h" +#include "podio_model/Measurement2DCollection.h" + + +namespace tdis::tracking { + + struct Measurement2DFactory : public JOmniFactory { + PodioInput m_mc_hits_in {this, {"TrackerHit"}}; + PodioOutput m_tracker_hits_out {this, "Measurement2D"}; + Service m_service_geometry{this}; + Parameter m_cfg_use_true_pos{this, "acts:use_true_position", true,"Use true hits xyz instead of digitized one"}; + + void Configure() { + m_service_geometry(); + } + + void ChangeRun(int32_t /*run_nr*/) { + } + + void Execute(int32_t /*run_nr*/, uint64_t /*evt_nr*/) { + using namespace Acts::UnitLiterals; + + + auto measurements = std::make_unique(); + auto plane_positions = m_service_geometry->GetPlanePositions(); + + + m_tracker_hits_out() = std::move(measurements); + } + }; +} \ No newline at end of file diff --git a/source/tdis/tracking/ReconstructedHitFactory.h b/source/tdis/tracking/ReconstructedHitFactory.h index 15cf2e5..42dedb5 100644 --- a/source/tdis/tracking/ReconstructedHitFactory.h +++ b/source/tdis/tracking/ReconstructedHitFactory.h @@ -5,6 +5,8 @@ #include "PadGeometryHelper.hpp" #include "podio_model/DigitizedMtpcMcHit.h" +#include "podio_model/Measurement2D.h" +#include "podio_model/Measurement2DCollection.h" #include "podio_model/MutableTrackerHit.h" #include "podio_model/TrackerHit.h" #include "podio_model/TrackerHitCollection.h" @@ -18,41 +20,50 @@ namespace { const double res = get_resolution(pixel_size); return res * res; } -} // namespace +} // namespace namespace tdis::tracking { struct ReconstructedHitFactory : public JOmniFactory { - PodioInput m_mc_hits_in {this, {"DigitizedMtpcMcHit"}}; - PodioOutput m_tracker_hits_out {this, "TrackerHit"}; + PodioInput m_mc_hits_in{this, {"DigitizedMtpcMcHit"}}; + PodioOutput m_tracker_hits_out{this, "TrackerHit"}; + PodioOutput m_measurements_out{this, "Measurement2D"}; Service m_service_geometry{this}; - Parameter m_cfg_use_true_pos{this, "acts:use_true_position", true,"Use true hits xyz instead of digitized one"}; + Service m_service_log{this}; + Parameter m_cfg_use_true_pos{this, "acts:use_true_position", true, + "Use true hits xyz instead of digitized one"}; + + std::shared_ptr m_log; void Configure() { m_service_geometry(); + m_log = m_service_log->logger("tracking:hit_reco"); } - void ChangeRun(int32_t /*run_nr*/) { - } + void ChangeRun(int32_t /*run_nr*/) {} - void Execute(int32_t /*run_nr*/, uint64_t /*evt_nr*/) { + void Execute(int32_t /*run_nr*/, uint64_t event_index) { using namespace Acts::UnitLiterals; - auto rec_hits = std::make_unique(); + auto measurements = std::make_unique(); auto plane_positions = m_service_geometry->GetPlanePositions(); - for (auto mc_hit : *m_mc_hits_in()) { + m_log->trace("ReconstructedHitFactory, reconstructing event: {}", event_index); + for (auto mc_hit : *m_mc_hits_in()) { const int plane = mc_hit.plane(); const int ring = mc_hit.ring(); const int pad = mc_hit.pad(); const double z_to_gem = mc_hit.zToGem(); - auto [pad_x,pad_y] = getPadCenter(ring, pad); - double plane_z = m_service_geometry().GetPlanePositions()[plane]; + m_log->trace("Plane {}, ring {}, pad {}, z_to_gem {}. True x {}, y {}, z {}", plane, ring, pad, z_to_gem, + mc_hit.truePosition().x, mc_hit.truePosition().y, mc_hit.truePosition().z); + auto [pad_x, pad_y] = getPadCenter(ring, pad); + double plane_z = plane_positions[plane]; + // Planes are set back to back like this: // | || || || | // This means that zToPlane is positive for even planes and negative for odd @@ -61,15 +72,13 @@ namespace tdis::tracking { // z = plane_z - zToPlane // for odd double calc_z = plane_z + (plane % 2 ? -z_to_gem : z_to_gem); - // Position edm4hep::Vector3f position; - if(m_cfg_use_true_pos() && !std::isnan(mc_hit.truePosition().x)) { + if (m_cfg_use_true_pos() && !std::isnan(mc_hit.truePosition().x)) { position.x = mc_hit.truePosition().x; position.y = mc_hit.truePosition().y; position.z = mc_hit.truePosition().z; - } - else { + } else { position.x = pad_x; position.y = pad_y; position.z = calc_z; @@ -80,18 +89,70 @@ namespace tdis::tracking { double xy_variance = get_variance(max_dimension); edm4eic::CovDiag3f cov{xy_variance, xy_variance, 1_cm}; // TODO correct Z variance - uint32_t cell_id = 1000000*mc_hit.plane() + 1000*mc_hit.ring() + mc_hit.pad(); + uint32_t cell_id = 1000000 * mc_hit.plane() + 1000 * mc_hit.ring() + mc_hit.pad(); + + auto hit + = rec_hits->create(cell_id, position, cov, static_cast(mc_hit.time()), + static_cast(1_ns), // TODO correct time resolution + static_cast(mc_hit.adc()), 0.0F); + + hit.rawHit(mc_hit); + + auto acts_det_element = m_service_geometry().GetDetectorElement(mc_hit.plane()); + + // Measurement 2d + + // variable surf_center not used anywhere; + + const auto& hit_pos = hit.position(); // 3d position + + Acts::Vector2 loc = Acts::Vector2::Zero(); + Acts::Vector2 pos; + auto onSurfaceTolerance + = 0.1 * Acts::UnitConstants::um; // By default, ACTS uses 0.1 micron as the on + // surface tolerance + + try { + // transform global position into local coordinates + // geometry context contains nothing here + pos = acts_det_element->surface() + .globalToLocal(Acts::GeometryContext(), + {hit_pos.x, hit_pos.y, plane_z}, + {0, 0, 0}, + onSurfaceTolerance) + .value(); + + loc[Acts::eBoundLoc0] = pos[0]; + loc[Acts::eBoundLoc1] = pos[1]; + } catch (std::exception& ex) { + m_log->warn( + "Can't convert globalToLocal for hit: plane={} ring={} pad={} RecoHit x={} y={} z={}. Reason: {}", + mc_hit.plane(), mc_hit.ring(), mc_hit.pad(), hit_pos.x, hit_pos.y, hit_pos.z, ex.what()); + continue; + } + + if (m_log->level() <= spdlog::level::trace) { + double surf_center_x = acts_det_element->surface().center(Acts::GeometryContext()).transpose()[0]; + double surf_center_y = acts_det_element->surface().center(Acts::GeometryContext()).transpose()[1]; + double surf_center_z = acts_det_element->surface().center(Acts::GeometryContext()).transpose()[2]; + m_log->trace(" hit position : {:>10.2f} {:>10.2f} {:>10.2f}", hit_pos.x, hit_pos.y, hit_pos.z); + m_log->trace(" surface center : {:>10.2f} {:>10.2f} {:>10.2f}", surf_center_x, surf_center_y, surf_center_z); + m_log->trace(" acts local center: {:>10.2f} {:>10.2f}", pos.transpose()[0], pos.transpose()[1]); + m_log->trace(" acts loc pos : {:>10.2f} {:>10.2f}", loc[Acts::eBoundLoc0], + loc[Acts::eBoundLoc1]); + } - rec_hits->create( - cell_id, - position, - cov, - static_cast(mc_hit.time()), - static_cast(1_ns), // TODO correct time resolution - static_cast(mc_hit.adc()), - 0.0F); + auto meas2D = measurements->create(); + meas2D.surface(acts_det_element->surface().geometryId().value()); // Surface for bound coordinates (geometryID) + meas2D.loc({static_cast(pos[0]), static_cast(pos[1])}); // 2D location on surface + meas2D.time(hit.time()); // Measurement time + // fixme: no off-diagonal terms. cov(0,1) = cov(1,0)?? + meas2D.covariance({cov(0, 0), cov(1, 1), hit.timeError() * hit.timeError(), cov(0, 1)}); // Covariance on location and time + meas2D.addweights(1.0); // Weight for each of the hits, mirrors hits array + meas2D.addhits(hit); } m_tracker_hits_out() = std::move(rec_hits); + m_measurements_out() = std::move(measurements); } }; -} \ No newline at end of file +} // namespace tdis::tracking \ No newline at end of file diff --git a/source/tdis/tracking/TruthTrackParameterFactory.h b/source/tdis/tracking/TruthTrackParameterFactory.h new file mode 100644 index 0000000..7cfa22f --- /dev/null +++ b/source/tdis/tracking/TruthTrackParameterFactory.h @@ -0,0 +1,97 @@ +#pragma once + +#include +#include + +#include "PadGeometryHelper.hpp" +#include "podio_model/DigitizedMtpcMcHit.h" +#include "podio_model/MutableTrackerHit.h" +#include "podio_model/TrackerHit.h" +#include "podio_model/TrackerHitCollection.h" + +namespace { + inline double get_resolution(const double pixel_size) { + constexpr const double sqrt_12 = 3.4641016151; + return pixel_size / sqrt_12; + } + inline double get_variance(const double pixel_size) { + const double res = get_resolution(pixel_size); + return res * res; + } +} // namespace + +namespace tdis::tracking { + + struct TruthTrackParameterFactory : public JOmniFactory { + PodioInput m_mc_hits_in {this, {"DigitizedMtpcMcHit"}}; + PodioOutput m_tracker_hits_out {this, "TrackerHit"}; + Service m_service_geometry{this}; + Parameter m_cfg_use_true_pos{this, "acts:use_true_position", true,"Use true hits xyz instead of digitized one"}; + + void Configure() { + m_service_geometry(); + } + + void ChangeRun(int32_t /*run_nr*/) { + } + + void Execute(int32_t /*run_nr*/, uint64_t /*evt_nr*/) { + using namespace Acts::UnitLiterals; + + + auto rec_hits = std::make_unique(); + auto plane_positions = m_service_geometry->GetPlanePositions(); + + for (auto mc_hit : *m_mc_hits_in()) { + + const int plane = mc_hit.plane(); + const int ring = mc_hit.ring(); + const int pad = mc_hit.pad(); + const double z_to_gem = mc_hit.zToGem(); + + auto [pad_x,pad_y] = getPadCenter(ring, pad); + double plane_z = m_service_geometry().GetPlanePositions()[plane]; + + + // Planes are set back to back like this: + // | || || || | + // This means that zToPlane is positive for even planes and negative for odd + // i.e. + // z = plane_z + zToPlane // for even + // z = plane_z - zToPlane // for odd + double calc_z = plane_z + (plane % 2 ? -z_to_gem : z_to_gem); + + + // Position + edm4hep::Vector3f position; + if(m_cfg_use_true_pos() && !std::isnan(mc_hit.truePosition().x)) { + position.x = mc_hit.truePosition().x; + position.y = mc_hit.truePosition().y; + position.z = mc_hit.truePosition().z; + } + else { + position.x = pad_x; + position.y = pad_y; + position.z = calc_z; + } + + // Covariance + double max_dimension = std::max(getPadApproxWidth(ring), getPadHight()); + double xy_variance = get_variance(max_dimension); + edm4eic::CovDiag3f cov{xy_variance, xy_variance, 1_cm}; // TODO correct Z variance + + uint32_t cell_id = 1000000*mc_hit.plane() + 1000*mc_hit.ring() + mc_hit.pad(); + + rec_hits->create( + cell_id, + position, + cov, + static_cast(mc_hit.time()), + static_cast(1_ns), // TODO correct time resolution + static_cast(mc_hit.adc()), + 0.0F); + } + m_tracker_hits_out() = std::move(rec_hits); + } + }; +} \ No newline at end of file From aee1808efa1bcf43c1873763f28ecdda169ddeba Mon Sep 17 00:00:00 2001 From: Dmitry Romanov Date: Tue, 19 Nov 2024 10:28:33 -0500 Subject: [PATCH 8/8] CKF Tracking --- source/tdis/CMakeLists.txt | 3 +- source/tdis/io/PodioWriteProcessor.hpp | 4 +- source/tdis/tdis_main.cpp | 10 +- source/tdis/tracking/CKFTracking.cc | 375 ++++++++++++++++++ source/tdis/tracking/CKFTracking.h | 91 +++++ .../tracking/TruthTrackParameterFactory.h | 163 +++++--- 6 files changed, 589 insertions(+), 57 deletions(-) create mode 100644 source/tdis/tracking/CKFTracking.cc create mode 100644 source/tdis/tracking/CKFTracking.h diff --git a/source/tdis/CMakeLists.txt b/source/tdis/CMakeLists.txt index 2c17dec..6b2323e 100644 --- a/source/tdis/CMakeLists.txt +++ b/source/tdis/CMakeLists.txt @@ -50,7 +50,8 @@ add_executable(tdis # tracking/Measurement2DFactory.h tracking/TruthTrackParameterFactory.h tracking/KalmanFittingFactory.h - + tracking/CKFTracking.h + tracking/CKFTracking.cc ) # ---------- FIND REQUIRED PACKAGES ------------- diff --git a/source/tdis/io/PodioWriteProcessor.hpp b/source/tdis/io/PodioWriteProcessor.hpp index 1278d51..ba11c75 100644 --- a/source/tdis/io/PodioWriteProcessor.hpp +++ b/source/tdis/io/PodioWriteProcessor.hpp @@ -75,7 +75,9 @@ class PodioWriteProcessor : public JEventProcessor { "DigitizedMtpcMcHit", // Digitized hits - "TrackerHit" + "TrackerHit", + "Measurement2D", + "TruthTrackInitParameters", }; PodioWriteProcessor(JApplication * app); diff --git a/source/tdis/tdis_main.cpp b/source/tdis/tdis_main.cpp index 0937994..33de276 100644 --- a/source/tdis/tdis_main.cpp +++ b/source/tdis/tdis_main.cpp @@ -9,9 +9,9 @@ // #include // #include +#include #include #include -#include #include @@ -21,6 +21,7 @@ #include "services/LogService.hpp" #include "tracking/ActsGeometryService.h" #include "tracking/ReconstructedHitFactory.h" +#include "tracking/TruthTrackParameterFactory.h" // #include "tracking/Measurement2DFactory.h" struct ProgramArguments { @@ -123,6 +124,13 @@ int main(int argc, char* argv[]) { {"TrackerHit", "Measurement2D"}); app.Add(reco_hit_generator); + auto truth_track_init_generator = new JOmniFactoryGeneratorT(); + truth_track_init_generator->AddWiring( + "TruthTrackParameterGenerator", + {"DigitizedMtpcMcTrack"}, + {"TruthTrackInitParameters"}); + app.Add(truth_track_init_generator); + // auto measurement_2d_generator = new JOmniFactoryGeneratorT(); // measurement_2d_generator->AddWiring("TrackerHitGenerator", {"TrackerHit"}, {"Measurement2D"}); // app.Add(measurement_2d_generator); diff --git a/source/tdis/tracking/CKFTracking.cc b/source/tdis/tracking/CKFTracking.cc new file mode 100644 index 0000000..93d33fe --- /dev/null +++ b/source/tdis/tracking/CKFTracking.cc @@ -0,0 +1,375 @@ +// SPDX-License-Identifier: LGPL-3.0-or-later +// Copyright (C) 2022 Whitney Armstrong, Wouter Deconinck, Dmitry Romanov, Shujie Li + +#include "CKFTracking.h" + +#include +#include +#include +#include +#if Acts_VERSION_MAJOR < 36 +#include +#endif +#include +#include +#if Acts_VERSION_MAJOR >= 32 +#include "Acts/EventData/ProxyAccessor.hpp" +#endif +#include +#include +#include +#include +#include +#include +#if Acts_VERSION_MAJOR >= 34 +#include "Acts/Propagator/AbortList.hpp" +#include "Acts/Propagator/EigenStepper.hpp" +#include "Acts/Propagator/MaterialInteractor.hpp" +#include "Acts/Propagator/Navigator.hpp" +#endif +#include +#if Acts_VERSION_MAJOR >= 34 +#include "Acts/Propagator/StandardAborters.hpp" +#endif +#include +#include +#include +#include +#include +#if Acts_VERSION_MAJOR >= 34 +#include "Acts/Utilities/TrackHelpers.hpp" +#endif +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include + +#include "ActsGeometryProvider.h" +#include "DD4hepBField.h" +#include "extensions/spdlog/SpdlogFormatters.h" // IWYU pragma: keep +#include "extensions/spdlog/SpdlogToActs.h" + +namespace eicrecon { + + using namespace Acts::UnitLiterals; + + // This array relates the Acts and EDM4eic covariance matrices, including + // the unit conversion to get from Acts units into EDM4eic units. + // + // Note: std::map is not constexpr, so we use a constexpr std::array + // std::array initialization need double braces since arrays are aggregates + // ref: https://en.cppreference.com/w/cpp/language/aggregate_initialization + static constexpr std::array, 6> edm4eic_indexed_units{{ + {Acts::eBoundLoc0, Acts::UnitConstants::mm}, + {Acts::eBoundLoc1, Acts::UnitConstants::mm}, + {Acts::eBoundPhi, 1.}, + {Acts::eBoundTheta, 1.}, + {Acts::eBoundQOverP, 1. / Acts::UnitConstants::GeV}, + {Acts::eBoundTime, Acts::UnitConstants::ns} + }}; + + CKFTracking::CKFTracking() { + } + + void CKFTracking::init(std::shared_ptr geo_svc, std::shared_ptr log) { + m_log = log; + m_acts_logger = eicrecon::getSpdlogLogger("CKF", m_log); + + m_geoSvc = geo_svc; + + m_BField = std::dynamic_pointer_cast(m_geoSvc->getFieldProvider()); + m_fieldctx = eicrecon::BField::BFieldVariant(m_BField); + + // eta bins, chi2 and #sourclinks per surface cutoffs + m_sourcelinkSelectorCfg = { + {Acts::GeometryIdentifier(), + {m_cfg.etaBins, m_cfg.chi2CutOff, + {m_cfg.numMeasurementsCutOff.begin(), m_cfg.numMeasurementsCutOff.end()} + } + }, + }; + m_trackFinderFunc = CKFTracking::makeCKFTrackingFunction(m_geoSvc->trackingGeometry(), m_BField, logger()); + } + + std::tuple< + std::vector, + std::vector + > + CKFTracking::process(const edm4eic::Measurement2DCollection& meas2Ds, + const edm4eic::TrackParametersCollection &init_trk_params) { + + + // create sourcelink and measurement containers + auto measurements = std::make_shared(); + + // need list here for stable addresses + std::list sourceLinkStorage; + ActsExamples::IndexSourceLinkContainer src_links; + src_links.reserve(meas2Ds.size()); + std::size_t hit_index = 0; + + + for (const auto& meas2D : meas2Ds) { + + // --follow example from ACTS to create source links + sourceLinkStorage.emplace_back(meas2D.getSurface(), hit_index); + ActsExamples::IndexSourceLink& sourceLink = sourceLinkStorage.back(); + // Add to output containers: + // index map and source link container are geometry-ordered. + // since the input is also geometry-ordered, new items can + // be added at the end. + src_links.insert(src_links.end(), sourceLink); + // --- + // Create ACTS measurements + Acts::Vector2 loc = Acts::Vector2::Zero(); + loc[Acts::eBoundLoc0] = meas2D.getLoc().a; + loc[Acts::eBoundLoc1] = meas2D.getLoc().b; + + + Acts::SquareMatrix2 cov = Acts::SquareMatrix2::Zero(); + cov(0, 0) = meas2D.getCovariance().xx; + cov(1, 1) = meas2D.getCovariance().yy; + cov(0, 1) = meas2D.getCovariance().xy; + cov(1, 0) = meas2D.getCovariance().xy; + +#if Acts_VERSION_MAJOR >= 36 + auto measurement = ActsExamples::makeFixedSizeMeasurement( + Acts::SourceLink{sourceLink}, loc, cov, Acts::eBoundLoc0, Acts::eBoundLoc1); +#else + auto measurement = Acts::makeMeasurement( + Acts::SourceLink{sourceLink}, loc, cov, Acts::eBoundLoc0, Acts::eBoundLoc1); +#endif + measurements->emplace_back(std::move(measurement)); + + hit_index++; + } + + ActsExamples::TrackParametersContainer acts_init_trk_params; + for (const auto& track_parameter: init_trk_params) { + + Acts::BoundVector params; + params(Acts::eBoundLoc0) = track_parameter.getLoc().a * Acts::UnitConstants::mm; // cylinder radius + params(Acts::eBoundLoc1) = track_parameter.getLoc().b * Acts::UnitConstants::mm; // cylinder length + params(Acts::eBoundPhi) = track_parameter.getPhi(); + params(Acts::eBoundTheta) = track_parameter.getTheta(); + params(Acts::eBoundQOverP) = track_parameter.getQOverP() / Acts::UnitConstants::GeV; + params(Acts::eBoundTime) = track_parameter.getTime() * Acts::UnitConstants::ns; + + double charge = std::copysign(1., track_parameter.getQOverP()); + + Acts::BoundSquareMatrix cov = Acts::BoundSquareMatrix::Zero(); + for (size_t i = 0; const auto& [a, x] : edm4eic_indexed_units) { + for (size_t j = 0; const auto& [b, y] : edm4eic_indexed_units) { + cov(a, b) = track_parameter.getCovariance()(i,j) * x * y; + ++j; + } + ++i; + } + + // Construct a perigee surface as the target surface + auto pSurface = Acts::Surface::makeShared(Acts::Vector3(0,0,0)); + + // Create parameters + acts_init_trk_params.emplace_back(pSurface, params, cov, Acts::ParticleHypothesis::pion()); + } + + //// Construct a perigee surface as the target surface + auto pSurface = Acts::Surface::makeShared(Acts::Vector3{0., 0., 0.}); + + ACTS_LOCAL_LOGGER(eicrecon::getSpdlogLogger("CKF", m_log, {"^No tracks found$"})); + + Acts::PropagatorPlainOptions pOptions; + pOptions.maxSteps = 10000; + + ActsExamples::PassThroughCalibrator pcalibrator; + ActsExamples::MeasurementCalibratorAdapter calibrator(pcalibrator, *measurements); + Acts::GainMatrixUpdater kfUpdater; +#if Acts_VERSION_MAJOR < 34 + Acts::GainMatrixSmoother kfSmoother; +#endif + Acts::MeasurementSelector measSel{m_sourcelinkSelectorCfg}; + + Acts::CombinatorialKalmanFilterExtensions + extensions; + extensions.calibrator.connect<&ActsExamples::MeasurementCalibratorAdapter::calibrate>( + &calibrator); + extensions.updater.connect< + &Acts::GainMatrixUpdater::operator()>( + &kfUpdater); +#if Acts_VERSION_MAJOR < 34 + extensions.smoother.connect< + &Acts::GainMatrixSmoother::operator()>( + &kfSmoother); +#endif + extensions.measurementSelector.connect< + &Acts::MeasurementSelector::select>( + &measSel); + + ActsExamples::IndexSourceLinkAccessor slAccessor; + slAccessor.container = &src_links; + Acts::SourceLinkAccessorDelegate + slAccessorDelegate; + slAccessorDelegate.connect<&ActsExamples::IndexSourceLinkAccessor::range>(&slAccessor); + + // Set the CombinatorialKalmanFilter options +#if Acts_VERSION_MAJOR < 34 + CKFTracking::TrackFinderOptions options( + m_geoctx, m_fieldctx, m_calibctx, slAccessorDelegate, + extensions, pOptions, &(*pSurface)); +#else + CKFTracking::TrackFinderOptions options( + m_geoctx, m_fieldctx, m_calibctx, slAccessorDelegate, + extensions, pOptions); +#endif + +#if Acts_VERSION_MAJOR >= 34 + Acts::Propagator, Acts::Navigator> extrapolator( + Acts::EigenStepper<>(m_BField), + Acts::Navigator({m_geoSvc->trackingGeometry()}, + logger().cloneWithSuffix("Navigator")), + logger().cloneWithSuffix("Propagator")); + + Acts::PropagatorOptions, + Acts::AbortList> + extrapolationOptions(m_geoctx, m_fieldctx); +#endif + + // Create track container + auto trackContainer = std::make_shared(); + auto trackStateContainer = std::make_shared(); + ActsExamples::TrackContainer acts_tracks(trackContainer, trackStateContainer); + + // Add seed number column + acts_tracks.addColumn("seed"); +#if Acts_VERSION_MAJOR >= 32 + Acts::ProxyAccessor seedNumber("seed"); +#else + Acts::TrackAccessor seedNumber("seed"); +#endif + + // Loop over seeds + for (std::size_t iseed = 0; iseed < acts_init_trk_params.size(); ++iseed) { + auto result = + (*m_trackFinderFunc)(acts_init_trk_params.at(iseed), options, acts_tracks); + + if (!result.ok()) { + m_log->debug("Track finding failed for seed {} with error {}", iseed, result.error()); + continue; + } + + // Set seed number for all found tracks + auto& tracksForSeed = result.value(); + for (auto& track : tracksForSeed) { + +#if Acts_VERSION_MAJOR >=34 + auto smoothingResult = Acts::smoothTrack(m_geoctx, track, logger()); + if (!smoothingResult.ok()) { + ACTS_ERROR("Smoothing for seed " + << iseed << " and track " << track.index() + << " failed with error " << smoothingResult.error()); + continue; + } + + auto extrapolationResult = Acts::extrapolateTrackToReferenceSurface( + track, *pSurface, extrapolator, extrapolationOptions, + Acts::TrackExtrapolationStrategy::firstOrLast, logger()); + if (!extrapolationResult.ok()) { + ACTS_ERROR("Extrapolation for seed " + << iseed << " and track " << track.index() + << " failed with error " << extrapolationResult.error()); + continue; + } +#endif + + seedNumber(track) = iseed; + } + } + + + // Move track states and track container to const containers + // NOTE Using the non-const containers leads to references to + // implicitly converted temporaries inside the Trajectories. + auto constTrackStateContainer = + std::make_shared( + std::move(*trackStateContainer)); + + auto constTrackContainer = + std::make_shared( + std::move(*trackContainer)); + + // FIXME JANA2 std::vector requires wrapping ConstTrackContainer, instead of: + //ConstTrackContainer constTracks(constTrackContainer, constTrackStateContainer); + std::vector constTracks_v; + constTracks_v.push_back( + new ActsExamples::ConstTrackContainer( + constTrackContainer, + constTrackStateContainer)); + auto& constTracks = *(constTracks_v.front()); + + // Seed number column accessor +#if Acts_VERSION_MAJOR >= 32 + const Acts::ConstProxyAccessor constSeedNumber("seed"); +#else + const Acts::ConstTrackAccessor constSeedNumber("seed"); +#endif + + // Prepare the output data with MultiTrajectory, per seed + std::vector acts_trajectories; + acts_trajectories.reserve(init_trk_params.size()); + + ActsExamples::Trajectories::IndexedParameters parameters; + std::vector tips; + + std::optional lastSeed; + for (const auto& track : constTracks) { + if (!lastSeed) { + lastSeed = constSeedNumber(track); + } + + if (constSeedNumber(track) != lastSeed.value()) { + // make copies and clear vectors + acts_trajectories.push_back(new ActsExamples::Trajectories( + constTracks.trackStateContainer(), + tips, parameters)); + + tips.clear(); + parameters.clear(); + } + + lastSeed = constSeedNumber(track); + + tips.push_back(track.tipIndex()); + parameters.emplace( + std::pair{track.tipIndex(), + ActsExamples::TrackParameters{track.referenceSurface().getSharedPtr(), + track.parameters(), track.covariance(), + track.particleHypothesis()}}); + } + + if (tips.empty()) { + m_log->info("Last trajectory is empty"); + } + + // last entry: move vectors + acts_trajectories.push_back(new ActsExamples::Trajectories( + constTracks.trackStateContainer(), + std::move(tips), std::move(parameters))); + + return std::make_tuple(std::move(acts_trajectories), std::move(constTracks_v)); + } + +} // namespace eicrecon diff --git a/source/tdis/tracking/CKFTracking.h b/source/tdis/tracking/CKFTracking.h new file mode 100644 index 0000000..df67aec --- /dev/null +++ b/source/tdis/tracking/CKFTracking.h @@ -0,0 +1,91 @@ +// SPDX-License-Identifier: LGPL-3.0-or-later +// Copyright (C) 2022 Whitney Armstrong, Wouter Deconinck + +#pragma once + +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include + +#include "CKFTrackingConfig.h" +#include "DD4hepBField.h" +#include "algorithms/interfaces/WithPodConfig.h" + +class ActsGeometryProvider; + +namespace eicrecon { + +/** Fitting algorithm implementation . + * + * \ingroup tracking + */ + + class CKFTracking: public WithPodConfig { + public: + + + /// Find function that takes the above parameters + /// @note This is separated into a virtual interface to keep compilation units + /// small + class CKFTrackingFunction { + public: + virtual ~CKFTrackingFunction() = default; + + virtual TrackFinderResult operator()(const ActsExamples::TrackParameters&, + const TrackFinderOptions&, + ActsExamples::TrackContainer&) const = 0; + }; + + /// Create the track finder function implementation. + /// The magnetic field is intentionally given by-value since the variantresults + /// contains shared_ptr anyways. + static std::shared_ptr makeCKFTrackingFunction( + std::shared_ptr trackingGeometry, + std::shared_ptr magneticField, + const Acts::Logger& logger); + + CKFTracking(); + + void init(std::shared_ptr geo_svc, std::shared_ptr log); + + std::tuple< + std::vector, + std::vector + > + process(const edm4eic::Measurement2DCollection& meas2Ds, + const edm4eic::TrackParametersCollection &init_trk_params); + + private: + std::shared_ptr m_log; + std::shared_ptr m_acts_logger{nullptr}; + std::shared_ptr m_trackFinderFunc; + std::shared_ptr m_geoSvc; + + std::shared_ptr m_BField = nullptr; + Acts::GeometryContext m_geoctx; + Acts::CalibrationContext m_calibctx; + Acts::MagneticFieldContext m_fieldctx; + + Acts::MeasurementSelector::Config m_sourcelinkSelectorCfg; + + /// Private access to the logging instance + const Acts::Logger& logger() const { return *m_acts_logger; } + }; + +} // namespace eicrecon::Reco diff --git a/source/tdis/tracking/TruthTrackParameterFactory.h b/source/tdis/tracking/TruthTrackParameterFactory.h index 7cfa22f..1840251 100644 --- a/source/tdis/tracking/TruthTrackParameterFactory.h +++ b/source/tdis/tracking/TruthTrackParameterFactory.h @@ -2,96 +2,151 @@ #include #include +#include +#include +#include +#include "ActsGeometryService.h" #include "PadGeometryHelper.hpp" #include "podio_model/DigitizedMtpcMcHit.h" +#include "podio_model/DigitizedMtpcMcTrack.h" #include "podio_model/MutableTrackerHit.h" +#include "podio_model/TrackParametersCollection.h" #include "podio_model/TrackerHit.h" #include "podio_model/TrackerHitCollection.h" -namespace { - inline double get_resolution(const double pixel_size) { - constexpr const double sqrt_12 = 3.4641016151; - return pixel_size / sqrt_12; - } - inline double get_variance(const double pixel_size) { - const double res = get_resolution(pixel_size); - return res * res; - } -} // namespace namespace tdis::tracking { struct TruthTrackParameterFactory : public JOmniFactory { - PodioInput m_mc_hits_in {this, {"DigitizedMtpcMcHit"}}; - PodioOutput m_tracker_hits_out {this, "TrackerHit"}; + PodioInput m_mc_tracks_in {this, {"DigitizedMtpcMcTrack"}}; + PodioOutput m_trackers_out {this, "TruthTrackInitParameters"}; Service m_service_geometry{this}; + Service m_service_log{this}; + Parameter m_cfg_use_true_pos{this, "acts:use_true_position", true,"Use true hits xyz instead of digitized one"}; + Parameter m_cfg_momentum_smear{this, "acts:track_init:momentum_smear", 0.1,"GeV, Momentum smear for truth track initialization"}; + std::shared_ptr m_log; + void Configure() { m_service_geometry(); + m_log = m_service_log->logger("tracking:hit_reco"); } void ChangeRun(int32_t /*run_nr*/) { } + // Function to generate the next value in a normal distribution + double generateNormal(double mean, double stddev) { + // Create a random device and a generator + static std::random_device rd; + static std::mt19937 generator(rd()); + + // Create a normal distribution with the given mean and standard deviation + std::normal_distribution distribution(mean, stddev); + + // Generate and return the next value + return distribution(generator); + } + void Execute(int32_t /*run_nr*/, uint64_t /*evt_nr*/) { using namespace Acts::UnitLiterals; - auto rec_hits = std::make_unique(); + auto track_parameters = std::make_unique(); auto plane_positions = m_service_geometry->GetPlanePositions(); - for (auto mc_hit : *m_mc_hits_in()) { - const int plane = mc_hit.plane(); - const int ring = mc_hit.ring(); - const int pad = mc_hit.pad(); - const double z_to_gem = mc_hit.zToGem(); + // Loop over input particles + for (const auto& mc_track: *m_mc_tracks_in()) { - auto [pad_x,pad_y] = getPadCenter(ring, pad); - double plane_z = m_service_geometry().GetPlanePositions()[plane]; + // We need at least one hit + if(!mc_track.hits_size()) { + continue; + } + // require close to interaction vertex + auto v = mc_track.hits().at(0).truePosition(); // Use it as a vertex for now + double vx = v.x; + double vy = v.y; + double vz = v.z; - // Planes are set back to back like this: - // | || || || | - // This means that zToPlane is positive for even planes and negative for odd - // i.e. - // z = plane_z + zToPlane // for even - // z = plane_z - zToPlane // for odd - double calc_z = plane_z + (plane % 2 ? -z_to_gem : z_to_gem); + double magnitude = mc_track.momentum(); + double theta = mc_track.theta(); + double phi = mc_track.phi(); + const auto eta = -std::log(std::tan(theta/2)); + double px = magnitude * std::sin(theta) * std::cos(phi); + double py = magnitude * std::sin(theta) * std::sin(phi); + double pz = magnitude * std::cos(theta); - // Position - edm4hep::Vector3f position; - if(m_cfg_use_true_pos() && !std::isnan(mc_hit.truePosition().x)) { - position.x = mc_hit.truePosition().x; - position.y = mc_hit.truePosition().y; - position.z = mc_hit.truePosition().z; - } - else { - position.x = pad_x; - position.y = pad_y; - position.z = calc_z; + // require minimum momentum + const auto& p = mc_track.momentum(); + const auto pmag = std::hypot(px, py, pz); + + + + // modify initial momentum to avoid bleeding truth to results when fit fails + const auto pinit = pmag * generateNormal(1, m_cfg_momentum_smear()*Acts::UnitConstants::GeV); + + // define line surface for local position values + auto perigee = Acts::Surface::makeShared(Acts::Vector3(0,0,0)); + + // track particle back to transverse point-of-closest approach + // with respect to the defined line surface + auto linesurface_parameter = -(vx*px + vy*py)/(px*px + py*py); + + auto xpca = v.x + linesurface_parameter*px; + auto ypca = v.y + linesurface_parameter*py; + auto zpca = v.z + linesurface_parameter*pz; + + Acts::Vector3 global(xpca, ypca, zpca); + Acts::Vector3 direction(sin(theta)*cos(phi), sin(theta)*sin(phi), cos(theta)); + + // convert from global to local coordinates using the defined line surface + auto local = perigee->globalToLocal(m_service_geometry->GetActsGeometryContext(), global, direction); + + if(!local.ok()) + { + m_log->error("skipping the track because globaltoLocal function failed"); + continue; } - // Covariance - double max_dimension = std::max(getPadApproxWidth(ring), getPadHight()); - double xy_variance = get_variance(max_dimension); - edm4eic::CovDiag3f cov{xy_variance, xy_variance, 1_cm}; // TODO correct Z variance - - uint32_t cell_id = 1000000*mc_hit.plane() + 1000*mc_hit.ring() + mc_hit.pad(); - - rec_hits->create( - cell_id, - position, - cov, - static_cast(mc_hit.time()), - static_cast(1_ns), // TODO correct time resolution - static_cast(mc_hit.adc()), - 0.0F); + Acts::Vector2 localpos = local.value(); + double charge = 1; // TODO ??? Proton? + + // Insert into edm4eic::TrackParameters, which uses numerical values in its specified units + auto track_parameter = track_parameters->create(); + track_parameter.type(-1); // type --> seed(-1) + track_parameter.loc({static_cast(localpos(0)), static_cast(localpos(1))}); // 2d location on surface [mm] + track_parameter.phi(phi); // phi [rad] + track_parameter.theta(theta); // theta [rad] + track_parameter.qOverP(charge / (pinit)); // Q/p [e/GeV] + track_parameter.time(mc_track.hits().at(0).time()); // time [ns] + edm4eic::Cov6f cov; + cov(0,0) = 1.0; // loc0 + cov(1,1) = 1.0; // loc1 + cov(2,2) = 0.05; // phi + cov(3,3) = 0.01; // theta + cov(4,4) = 0.1; // qOverP + cov(5,5) = 10e9; // time + track_parameter.covariance(cov); + + // // Debug output + // if (m_log->level() <= spdlog::level::debug) { + // m_log->debug("Invoke track finding seeded by truth particle with:"); + // m_log->debug(" p = {} GeV (smeared to {} GeV)", pmag / dd4hep::GeV, pinit / dd4hep::GeV); + // m_log->debug(" q = {}", charge); + // m_log->debug(" q/p = {} e/GeV (smeared to {} e/GeV)", charge / (pmag / dd4hep::GeV), charge / (pinit / dd4hep::GeV)); + // m_log->debug(" theta = {}", theta); + // m_log->debug(" phi = {}", phi); + // } } - m_tracker_hits_out() = std::move(rec_hits); + + + + m_trackers_out() = std::move(track_parameters); } }; } \ No newline at end of file