algorithms/fardetectors/FarDetectorLinearProjection.cc

0001 // SPDX-License-Identifier: LGPL-3.0-or-later
0002 // Copyright (C) 2023-2025, Simon Gardner
0003
0004 #include <edm4eic/Cov6f.h>
0005 #include <edm4eic/TrackCollection.h>
0006 #include <edm4eic/vector_utils.h>
0007 #include <edm4hep/Vector2f.h>
0008 #include <edm4hep/Vector3f.h>
0009 #include <fmt/core.h>
0010 #include <Eigen/LU>
0011 #include <cstdint>
0012 #include <gsl/pointers>
0013 #include <vector>
0014
0015 #include "algorithms/fardetectors/FarDetectorLinearProjection.h"
0016 #include "algorithms/fardetectors/FarDetectorLinearProjectionConfig.h"
0017
0018 namespace eicrecon {
0019
0020 void FarDetectorLinearProjection::init() {
0021
0022   // plane position
0023   m_plane_position << m_cfg.plane_position[0], m_cfg.plane_position[1], m_cfg.plane_position[2];
0024   m_directions.block<3, 1>(0, 0) << m_cfg.plane_a[0], m_cfg.plane_a[1], m_cfg.plane_a[2];
0025   m_directions.block<3, 1>(0, 1) << m_cfg.plane_b[0], m_cfg.plane_b[1], m_cfg.plane_b[2];
0026 }
0027
0028 void FarDetectorLinearProjection::process(const FarDetectorLinearProjection::Input& input,
0029                                           const FarDetectorLinearProjection::Output& output) const {
0030
0031   const auto [inputTracks] = input;
0032   auto [outputTracks]      = output;
0033
0034   Eigen::Matrix3d directions = m_directions;
0035
0036   for (const auto& track : *inputTracks) {
0037
0038     Eigen::Vector3d point_position(track.getPosition().x, track.getPosition().y,
0039                                    track.getPosition().z);
0040     Eigen::Vector3d positionDiff = point_position - m_plane_position;
0041
0042     // Convert spherical coordinates to Cartesian
0043     double x = track.getMomentum().x;
0044     double y = track.getMomentum().y;
0045     double z = track.getMomentum().z;
0046     directions.block<3, 1>(0, 2) << x, y, z;
0047
0048     auto projectedPoint = directions.inverse() * positionDiff;
0049
0050     // Create track parameters edm4eic structure
0051     // TODO - populate more of the fields correctly
0052     std::int32_t type = 0;
0053     // Surface ID not used in this context
0054     std::uint64_t surface = 0;
0055     // Plane Point
0056     edm4hep::Vector2f loc(projectedPoint[0], projectedPoint[1]); //Temp unit transform
0057     float theta     = edm4eic::anglePolar(track.getMomentum());
0058     float phi       = edm4eic::angleAzimuthal(track.getMomentum());
0059     float qOverP    = 0.;
0060     float time      = 0;
0061     int32_t pdgCode = 11;
0062     // Point Error
0063     edm4eic::Cov6f error;
0064
0065     debug("Position:      a={},   b={}", loc.a, loc.b);
0066     debug("Direction: theta={}, phi={}", theta, phi);
0067
0068     outputTracks->create(type, surface, loc, theta, phi, qOverP, time, pdgCode, error);
0069   }
0070 }
0071
0072 } // namespace eicrecon