diff --git a/.codespell-ignore b/.codespell-ignore index 3493f94755..f27b1d0e08 100644 --- a/.codespell-ignore +++ b/.codespell-ignore @@ -2,3 +2,5 @@ gaus mapp te ths +millepede +mille diff --git a/src/algorithms/tracking/CMakeLists.txt b/src/algorithms/tracking/CMakeLists.txt index 80d4927583..e60b40dc00 100644 --- a/src/algorithms/tracking/CMakeLists.txt +++ b/src/algorithms/tracking/CMakeLists.txt @@ -27,3 +27,6 @@ plugin_include_directories(${PLUGIN_NAME} PUBLIC # Add libraries (same as target_include_directories but for both plugin and # library) plugin_link_libraries(${PLUGIN_NAME} Eigen3::Eigen particle_service_library) + +# MeasurementToMille.cc uses edm4eic::AlignmentDerivativeSet, available in +# edm4eic >= 8.10.0. diff --git a/src/algorithms/tracking/MeasurementToMille.cc b/src/algorithms/tracking/MeasurementToMille.cc new file mode 100644 index 0000000000..19b9220796 --- /dev/null +++ b/src/algorithms/tracking/MeasurementToMille.cc @@ -0,0 +1,157 @@ +// SPDX-License-Identifier: LGPL-3.0-or-later +// Copyright (C) 2025 ePIC Collaboration + +#include +#include +#include +#include +#include +#include +#include +#include +#include + +#include "MeasurementToMille.h" +#include "SiliconAlignmentLabels.h" +#include "algorithms/interfaces/ActsSvc.h" + +namespace eicrecon { + +void MeasurementToMille::init() { + m_geo = algorithms::ActsSvc::instance().acts_geometry_provider(); + + const auto& trackingGeo = *m_geo->trackingGeometry(); + const auto& dd4hepDetector = *m_geo->dd4hepDetector(); + + // Pass nullptr for the optional spdlog logger; the algorithms framework + // logger (info/debug methods below) is used instead. + m_surfaceToLayer = buildSiliconSurfaceLabelMap(trackingGeo, dd4hepDetector); + + info("Initialized: {} silicon surfaces mapped across {} known layers", m_surfaceToLayer.size(), + kNSiliconLayers); +} + +void MeasurementToMille::process(const Input& input, const Output& output) const { + const auto [tracks, measurements] = input; + auto [derivatives] = output; + + std::size_t nDerivatives = 0; + + for (const auto& track : *tracks) { + // ------------------------------------------------------------------ + // Track quality cuts + // ------------------------------------------------------------------ + const std::uint32_t ndf = track.getNdf(); + if (ndf == 0) { + continue; + } + const float chi2PerNDF = track.getChi2() / static_cast(ndf); + if (chi2PerNDF > m_cfg.maxChi2PerNDF) { + debug("Skipping track: chi2/NDF={:.2f} > {:.2f}", chi2PerNDF, m_cfg.maxChi2PerNDF); + continue; + } + + const auto& mom = track.getMomentum(); + const float p = edm4hep::utils::magnitude(mom); + if (p < m_cfg.minMomentum) { + debug("Skipping track: p={:.3f} GeV/c < {:.3f} GeV/c", p, m_cfg.minMomentum); + continue; + } + + // ------------------------------------------------------------------ + // Iterate over measurements associated with this track + // ------------------------------------------------------------------ + for (const auto& meas : track.getMeasurements()) { + const std::uint64_t geoId = meas.getSurface(); + + // Look up which silicon layer this surface belongs to + auto layerIt = m_surfaceToLayer.find(geoId); + if (layerIt == m_surfaceToLayer.end()) { + continue; // not a silicon alignment surface + } + const int layerIndex = layerIt->second; + + // Skip fixed (reference) layers + bool isFixed = false; + for (const int fl : m_cfg.fixedLayers) { + if (fl == layerIndex) { + isFixed = true; + break; + } + } + if (isFixed) { + continue; + } + + // ------------------------------------------------------------------ + // Residual and uncertainty + // + // TODO: Replace with the proper residual from the ACTS Kalman smoother: + // residual = measured_local_u - predicted_local_u + // where predicted_local_u comes from Acts::detail::makeTrackAlignmentState(). + // The current approximation sets predicted = 0, so residual = measured_u. + // This is only valid for perfectly centred sensors and is a placeholder + // until the ACTS Alignment kernel is accessible from EICrecon. + // ------------------------------------------------------------------ + const float measuredU = meas.getLoc().a; + const float residual = measuredU; // placeholder: predicted = 0 + const float residualUncert = std::sqrt(std::abs(meas.getCovariance().xx)); + + // ------------------------------------------------------------------ + // Local derivatives: ∂residual/∂local_track_parameters + // + // For a 1-D strip-like measurement in the bending plane (local u), + // the simplified Jacobian row is [1, 0, 0, 0, 0] corresponding to the + // five bound track parameters (loc0, loc1, phi, theta, qOverP). + // + // TODO: Use Acts::detail::makeTrackAlignmentState().localDerivatives + // for the correct Jacobian once ACTS integration is complete. + // ------------------------------------------------------------------ + constexpr int kNLocalParams = 5; + const std::array localDeriv = {1.f, 0.f, 0.f, 0.f, 0.f}; + + // ------------------------------------------------------------------ + // Global derivatives: ∂residual/∂alignment_DOFs + // + // For a translation in the sensor local-u direction (tx), moving the + // sensor by δ shifts the residual by −δ → derivative = -1. + // All other DOFs produce zero to first order at normal incidence. + // + // Layout: [tx, ty, tz, rx, ry, rz] + // + // TODO: Replace with Acts::detail::makeTrackAlignmentState().globalDerivatives + // which accounts for sensor orientation and track incidence angle. + // ------------------------------------------------------------------ + const std::array globalDeriv = { + -1.f, // tx: translation along local u directly shifts residual + 0.f, // ty: orthogonal to measurement direction + 0.f, // tz: along surface normal — negligible at normal incidence + 0.f, // rx: rotation about local x — no first-order effect on u at normal incidence + 0.f, // ry: rotation about local y — no first-order effect on u at normal incidence + 0.f, // rz: rotation about local z — couples v into u; zero in simplified model + }; + + // ------------------------------------------------------------------ + // Fill output object + // ------------------------------------------------------------------ + auto deriv = derivatives->create(geoId, residual, residualUncert); + + for (int i = 0; i < kNLocalParams; ++i) { + deriv.addToLocalDerivatives(localDeriv[i]); + } + + for (int dof = 0; dof < kAlignNDOF; ++dof) { + deriv.addToGlobalLabels(siliconAlignmentLabel(layerIndex, static_cast(dof))); + deriv.addToGlobalDerivatives(globalDeriv[dof]); + } + + ++nDerivatives; + } // end loop over track measurements + } // end loop over tracks + + debug("Filled {} AlignmentDerivativeSet entries from {} tracks", nDerivatives, tracks->size()); + // The measurements collection is available for future use (e.g. outlier cross-checks) + (void)measurements; +} + +} // namespace eicrecon diff --git a/src/algorithms/tracking/MeasurementToMille.h b/src/algorithms/tracking/MeasurementToMille.h new file mode 100644 index 0000000000..50d98abcfd --- /dev/null +++ b/src/algorithms/tracking/MeasurementToMille.h @@ -0,0 +1,59 @@ +// SPDX-License-Identifier: LGPL-3.0-or-later +// Copyright (C) 2025 ePIC Collaboration + +#pragma once + +#include +#include +#include +#include +#include +#include +#include +#include +#include + +#include "ActsGeometryProvider.h" +#include "MeasurementToMilleConfig.h" +#include "algorithms/interfaces/WithPodConfig.h" + +namespace eicrecon { + +using MeasurementToMilleAlgorithm = algorithms::Algorithm< + algorithms::Input, + algorithms::Output>; + +/// Converts reconstructed tracks and measurements into Millepede-II alignment +/// derivative sets for silicon tracker alignment. +/// +/// For each track passing quality cuts the algorithm iterates over the +/// Measurement2D hits associated with the track, looks up the alignment layer +/// from the ACTS geometry, and fills one AlignmentDerivativeSet per +/// (track, surface) pair with approximate local and global derivatives. +/// +/// NOTE: The current implementation uses simplified (approximate) derivatives. +/// The correct derivatives require calling Acts::detail::makeTrackAlignmentState() +/// on the ACTS track states. Once the ACTS Alignment kernel is accessible from +/// EICrecon, replace the simplified code with the proper computation. +class MeasurementToMille : public MeasurementToMilleAlgorithm, + public WithPodConfig { +public: + MeasurementToMille(std::string_view name) + : MeasurementToMilleAlgorithm{ + name, + {"inputTracks", "inputMeasurements"}, + {"outputAlignmentDerivatives"}, + "Fills Millepede-II alignment derivatives for silicon tracker"} {} + + void init() final; + void process(const Input&, const Output&) const final; + +private: + /// Map from Acts::GeometryIdentifier (as uint64_t) to 0-based silicon layer index. + /// Built once during init() by walking the ACTS TrackingGeometry. + std::unordered_map m_surfaceToLayer; + + std::shared_ptr m_geo; +}; + +} // namespace eicrecon diff --git a/src/algorithms/tracking/MeasurementToMilleConfig.h b/src/algorithms/tracking/MeasurementToMilleConfig.h new file mode 100644 index 0000000000..9e25ffca50 --- /dev/null +++ b/src/algorithms/tracking/MeasurementToMilleConfig.h @@ -0,0 +1,21 @@ +// SPDX-License-Identifier: LGPL-3.0-or-later +// Copyright (C) 2025 ePIC Collaboration + +#pragma once + +#include + +namespace eicrecon { + +/// Configuration for MeasurementToMille: selects tracks suitable for +/// silicon tracker alignment and controls which layers are constrained. +struct MeasurementToMilleConfig { + /// Maximum chi2/NDF for a track to be used in alignment + float maxChi2PerNDF = 5.0f; + /// Minimum track momentum [GeV/c] + float minMomentum = 1.0f; + /// 0-based layer indices to fix (skip) in Millepede, e.g. {0} for a reference layer + std::vector fixedLayers; +}; + +} // namespace eicrecon diff --git a/src/algorithms/tracking/SiliconAlignmentLabels.h b/src/algorithms/tracking/SiliconAlignmentLabels.h new file mode 100644 index 0000000000..fa3b56cd07 --- /dev/null +++ b/src/algorithms/tracking/SiliconAlignmentLabels.h @@ -0,0 +1,177 @@ +// SPDX-License-Identifier: LGPL-3.0-or-later +// Copyright (C) 2025 ePIC Collaboration + +#pragma once + +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include + +// Version-dependent ACTS DD4hep include +#if __has_include() +#include +using DD4hepDetElem = ActsPlugins::DD4hepDetectorElement; +#else +#include +using DD4hepDetElem = Acts::DD4hepDetectorElement; +#endif + +namespace eicrecon { + +/// Millepede alignment degrees of freedom (per layer, local sensor frame). +enum class AlignmentDOF : int { + kTx = 0, ///< Translation along local x [mm] + kTy = 1, ///< Translation along local y [mm] + kTz = 2, ///< Translation along local z [mm] + kRx = 3, ///< Rotation about local x [rad] + kRy = 4, ///< Rotation about local y [rad] + kRz = 5, ///< Rotation about local z [rad] + kNDOF = 6, +}; + +/// Number of DOFs per alignable layer. +static constexpr int kAlignNDOF = static_cast(AlignmentDOF::kNDOF); + +/// Ordered list of silicon tracker layer path prefixes (DD4hep DetElement paths). +/// +/// The order defines the Millepede label assignment: +/// label = layer_index * kAlignNDOF + dof_index + 1 (1-based positive integer) +/// +/// Layer ordering matches calibrations/alignment/silicon_misalignment.xml. +/// Each entry is the *prefix* of the DetElement path so that individual +/// sensitive elements (modules, sensors) deeper in the hierarchy are matched +/// to their parent layer. +static constexpr std::array kSiliconLayerPaths = {{ + // Layer 0: labels 1 - 6 + "/world/MiddleSiBarrelSubAssembly/SagittaSiBarrel/SagittaSiBarrel_layer1", + // Layer 1: labels 7 - 12 + "/world/OuterSiBarrelSubAssembly/OuterSiBarrel/OuterSiBarrel_layer2", + // Layer 2: labels 13 - 18 + "/world/InnerSiTrackerSubAssembly/InnerTrackerEndcapP/InnerTrackerEndcapP_layer1_P", + // Layer 3: labels 19 - 24 + "/world/InnerSiTrackerSubAssembly/InnerTrackerEndcapN/InnerTrackerEndcapN_layer1_N", + // Layer 4: labels 25 - 30 + "/world/MiddleSiEndcapSubAssembly/MiddleTrackerEndcapP/MiddleTrackerEndcapP_layer1_P", + // Layer 5: labels 31 - 36 + "/world/MiddleSiEndcapSubAssembly/MiddleTrackerEndcapN/MiddleTrackerEndcapN_layer1_N", + // Layer 6: labels 37 - 42 + "/world/OuterSiEndcapSubAssembly/OuterTrackerEndcapP/OuterTrackerEndcapP_layer2_P", + // Layer 7: labels 43 - 48 + "/world/OuterSiEndcapSubAssembly/OuterTrackerEndcapP/OuterTrackerEndcapP_layer3_P", + // Layer 8: labels 49 - 54 + "/world/OuterSiEndcapSubAssembly/OuterTrackerEndcapP/OuterTrackerEndcapP_layer4_P", + // Layer 9: labels 55 - 60 + "/world/OuterSiEndcapSubAssembly/OuterTrackerEndcapN/OuterTrackerEndcapN_layer2_N", + // Layer 10: labels 61 - 66 + "/world/OuterSiEndcapSubAssembly/OuterTrackerEndcapN/OuterTrackerEndcapN_layer3_N", + // Layer 11: labels 67 - 72 + "/world/OuterSiEndcapSubAssembly/OuterTrackerEndcapN/OuterTrackerEndcapN_layer4_N", +}}; + +/// Human-readable name for each layer (same order as kSiliconLayerPaths). +static constexpr std::array kSiliconLayerNames = {{ + "SagittaSiBarrel_layer1", + "OuterSiBarrel_layer2", + "InnerTrackerEndcapP_layer1_P", + "InnerTrackerEndcapN_layer1_N", + "MiddleTrackerEndcapP_layer1_P", + "MiddleTrackerEndcapN_layer1_N", + "OuterTrackerEndcapP_layer2_P", + "OuterTrackerEndcapP_layer3_P", + "OuterTrackerEndcapP_layer4_P", + "OuterTrackerEndcapN_layer2_N", + "OuterTrackerEndcapN_layer3_N", + "OuterTrackerEndcapN_layer4_N", +}}; + +static constexpr std::size_t kNSiliconLayers = kSiliconLayerPaths.size(); + +/// Return the Millepede base label (DOF 0 label) for a given 0-based layer index. +/// All DOF labels for this layer are [baseLabel, baseLabel + kAlignNDOF - 1]. +inline int siliconAlignmentBaseLabel(int layerIndex) { return layerIndex * kAlignNDOF + 1; } + +/// Return the Millepede label for a specific (layer, DOF) pair. +inline int siliconAlignmentLabel(int layerIndex, AlignmentDOF dof) { + return siliconAlignmentBaseLabel(layerIndex) + static_cast(dof); +} + +/// Decode a Millepede label back to (layerIndex, dof). +/// Returns std::nullopt if the label is out of range. +inline std::optional> decodeSiliconAlignmentLabel(int label) { + if (label < 1 || label > static_cast(kNSiliconLayers) * kAlignNDOF) { + return std::nullopt; + } + int zero_based = label - 1; + int layerIndex = zero_based / kAlignNDOF; + auto dof = static_cast(zero_based % kAlignNDOF); + return std::make_pair(layerIndex, dof); +} + +/// Build a map from Acts::GeometryIdentifier (encoded as uint64_t) to the +/// 0-based silicon layer index by walking the ACTS TrackingGeometry and +/// matching each sensitive surface's DD4hep DetElement path against the +/// kSiliconLayerPaths prefixes. +/// +/// The map is built once at algorithm initialization and queried per surface. +/// Returns only surfaces belonging to the 12 known silicon alignment layers. +inline std::unordered_map +buildSiliconSurfaceLabelMap(const Acts::TrackingGeometry& trackingGeo, + const dd4hep::Detector& dd4hepDetector, + std::shared_ptr log = nullptr) { + std::unordered_map surfaceToLayer; + + auto volman = dd4hepDetector.volumeManager(); + + trackingGeo.visitSurfaces([&](const Acts::Surface* surface) { + if (surface == nullptr) { + return; + } + +#if Acts_VERSION_MAJOR >= 45 + const auto* det_element = dynamic_cast(surface->surfacePlacement()); +#else + const auto* det_element = + dynamic_cast(surface->associatedDetectorElement()); +#endif + if (det_element == nullptr) { + return; + } + + auto* vol_ctx = volman.lookupContext(det_element->identifier()); + if (vol_ctx == nullptr) { + return; + } + const std::string path = vol_ctx->element.path(); + + // Find which silicon alignment layer this surface belongs to by + // checking whether the path starts with one of the known layer prefixes. + for (std::size_t i = 0; i < kNSiliconLayers; ++i) { + if (path.rfind(kSiliconLayerPaths[i], 0) == 0) { + const auto geoId = surface->geometryId().value(); + surfaceToLayer.emplace(geoId, static_cast(i)); + if (log) { + log->debug("SiliconAlignmentLabels: surface 0x{:016x} → layer {} ({}) path: {}", geoId, i, + kSiliconLayerNames[i], path); + } + break; + } + } + }); + + if (log) { + log->info("SiliconAlignmentLabels: mapped {} surfaces to {} layers", surfaceToLayer.size(), + kNSiliconLayers); + } + return surfaceToLayer; +} + +} // namespace eicrecon diff --git a/src/factories/tracking/MeasurementToMille_factory.h b/src/factories/tracking/MeasurementToMille_factory.h new file mode 100644 index 0000000000..2bf48a4f1c --- /dev/null +++ b/src/factories/tracking/MeasurementToMille_factory.h @@ -0,0 +1,54 @@ +// SPDX-License-Identifier: LGPL-3.0-or-later +// Copyright (C) 2025 ePIC Collaboration + +#pragma once + +#include +#include +#include +#include +#include +#include +#include +#include + +#include "algorithms/tracking/MeasurementToMille.h" +#include "algorithms/tracking/MeasurementToMilleConfig.h" +#include "extensions/jana/JOmniFactory.h" + +namespace eicrecon { + +class MeasurementToMille_factory + : public JOmniFactory { + +private: + using AlgoT = eicrecon::MeasurementToMille; + std::unique_ptr m_algo; + + PodioInput m_tracks_input{this}; + PodioInput m_measurements_input{this}; + PodioOutput m_derivatives_output{this}; + + ParameterRef m_maxChi2PerNDF{this, "maxChi2PerNDF", config().maxChi2PerNDF, + "Maximum chi2/NDF for tracks used in alignment"}; + ParameterRef m_minMomentum{this, "minMomentum", config().minMomentum, + "Minimum track momentum [GeV/c]"}; + ParameterRef> m_fixedLayers{this, "fixedLayers", config().fixedLayers, + "0-based layer indices to fix in Millepede"}; + + Service m_algorithmsInit{this}; + +public: + void Configure() { + m_algo = std::make_unique(this->GetPrefix()); + m_algo->level(static_cast(logger()->level())); + m_algo->applyConfig(config()); + m_algo->init(); + } + + void Process(int32_t /* run_number */, uint64_t /* event_number */) { + m_algo->process({m_tracks_input(), m_measurements_input()}, {m_derivatives_output().get()}); + } +}; + +} // namespace eicrecon diff --git a/src/global/CMakeLists.txt b/src/global/CMakeLists.txt index abe32a756d..5f1b4fe73e 100644 --- a/src/global/CMakeLists.txt +++ b/src/global/CMakeLists.txt @@ -1,4 +1,5 @@ add_subdirectory(tracking) +add_subdirectory(alignment) add_subdirectory(reco) add_subdirectory(pid) add_subdirectory(pid_lut) diff --git a/src/global/alignment/CMakeLists.txt b/src/global/alignment/CMakeLists.txt new file mode 100644 index 0000000000..cd80dd654d --- /dev/null +++ b/src/global/alignment/CMakeLists.txt @@ -0,0 +1,25 @@ +get_filename_component(PLUGIN_NAME ${CMAKE_CURRENT_LIST_DIR} NAME) + +# Creates ${PLUGIN_NAME}_plugin and ${PLUGIN_NAME}_library targets with default +# includes, libraries, and installation paths. +plugin_add(${PLUGIN_NAME} PLUGIN_USE_CC_ONLY) + +# Pick up alignment.cc (and any future sources) automatically. +plugin_glob_all(${PLUGIN_NAME}) + +# Dependencies shared with the tracking plugin +plugin_add_dd4hep(${PLUGIN_NAME}) +plugin_add_acts(${PLUGIN_NAME}) +plugin_add_event_model(${PLUGIN_NAME}) + +# Link against the tracking algorithm library so MeasurementToMille is +# available. +plugin_link_libraries(${PLUGIN_NAME} algorithms_tracking_library) + +# NOTE: The edm4eic::AlignmentDerivativeSet type is defined in the alignment +# branch of EDM4eic (~/git/EDM4eic/build-align/edm4eic/). Once it is merged +# into an official EDM4eic release and found by find_package(EDM4EIC), no extra +# include path is needed. Until then, add the build directory manually: +# +# plugin_include_directories(${PLUGIN_NAME} SYSTEM PUBLIC +# $ENV{HOME}/git/EDM4eic/build-align) diff --git a/src/global/alignment/alignment.cc b/src/global/alignment/alignment.cc new file mode 100644 index 0000000000..92af0aa8bd --- /dev/null +++ b/src/global/alignment/alignment.cc @@ -0,0 +1,50 @@ +// SPDX-License-Identifier: LGPL-3.0-or-later +// Copyright (C) 2025 ePIC Collaboration + +#include +#include +#include +#include +#include +#include +#include + +#include "extensions/jana/JOmniFactoryGeneratorT.h" +#include "factories/tracking/MeasurementToMille_factory.h" + +extern "C" { +void InitPlugin(JApplication* app) { + InitJANAPlugin(app); + + using namespace eicrecon; + + // Register the MeasurementToMille factory. + // + // Input collections: + // - "CentralCKFTracks" (edm4eic::TrackCollection) + // - "CentralTrackerMeasurements" (edm4eic::Measurement2DCollection) + // + // Output collection: + // - "SiliconAlignmentDerivatives" (edm4eic::AlignmentDerivativeSetCollection) + // + // NOTE: The EDM4eic AlignmentDerivativeSet type used here is defined in the + // ePIC alignment branch of EDM4eic and is not yet part of an official EDM4eic + // release. Update the EDM4eic dependency to the alignment-enabled build before + // compiling this plugin. + app->Add(new JOmniFactoryGeneratorT( + "SiliconAlignmentDerivatives", + { + "CentralCKFTracks", + "CentralTrackerMeasurements", + }, + { + "SiliconAlignmentDerivatives", + }, + { + .maxChi2PerNDF = 5.0f, + .minMomentum = 1.0f, + .fixedLayers = {0}, // fix SagittaSiBarrel_layer1 as reference by default + }, + app)); +} +} // extern "C"