algorithms/fardetectors/FarDetectorLinearProjection.cc

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