algorithms/fardetectors/FarDetectorLinearTracking.cc

0001 // SPDX-License-Identifier: LGPL-3.0-or-later
0002 // Copyright (C) 2023, Simon Gardner
0003
0004 #include <edm4eic/Cov2f.h>
0005 #include <edm4eic/Cov3f.h>
0006 #include <edm4eic/TrackPoint.h>
0007 #include <edm4eic/TrackSegmentCollection.h>
0008 #include <edm4eic/vector_utils.h>
0009 #include <edm4hep/TrackerHitCollection.h>
0010 #include <edm4hep/Vector3d.h>
0011 #include <edm4hep/Vector3f.h>
0012 #include <edm4hep/utils/vector_utils.h>
0013 #include <fmt/core.h>
0014 #include <stdint.h>
0015 #include <Eigen/Geometry>
0016 #include <Eigen/Householder>
0017 #include <Eigen/Jacobi>
0018 #include <Eigen/QR>
0019 #include <Eigen/SVD>
0020 #include <cmath>
0021 #include <utility>
0022
0023 #include "FarDetectorLinearTracking.h"
0024 #include "algorithms/fardetectors/FarDetectorLinearTrackingConfig.h"
0025
0026 namespace eicrecon {
0027
0028     void FarDetectorLinearTracking::init() {
0029
0030       // For changing how strongly each layer hit is in contributing to the fit
0031       m_layerWeights = Eigen::VectorXd::Constant(m_cfg.n_layer,1);
0032
0033       // For checking the direction of the track from theta and phi angles
0034       m_optimumDirection = Eigen::Vector3d::UnitZ();
0035       m_optimumDirection = Eigen::AngleAxisd(m_cfg.optimum_theta,Eigen::Vector3d::UnitY())*m_optimumDirection;
0036       m_optimumDirection = Eigen::AngleAxisd(m_cfg.optimum_phi,Eigen::Vector3d::UnitZ())*m_optimumDirection;
0037
0038     }
0039
0040     void FarDetectorLinearTracking::process(
0041         const FarDetectorLinearTracking::Input& input,
0042         const FarDetectorLinearTracking::Output& output) const {
0043
0044         const auto [inputhits] = input;
0045         auto [outputTracks] = output;
0046
0047         // Check the number of input collections is correct
0048         int nCollections = inputhits.size();
0049         if(nCollections!=m_cfg.n_layer){
0050           error("Wrong number of input collections passed to algorithm");
0051           return;
0052         }
0053
0054         // Check there aren't too many hits in any layer to handle
0055         // Temporary limit of number of hits per layer before Kalman filtering/GNN implemented
0056         // TODO - Implement more sensible solution
0057         for(const auto& layerHits: inputhits){
0058           if((*layerHits).size()>m_cfg.layer_hits_max){
0059             info("Too many hits in layer");
0060             return;
0061           }
0062         }
0063
0064         // Create a matrix to store the hit positions
0065         Eigen::MatrixXd hitMatrix(3,m_cfg.n_layer);
0066         // Loop over all combinations of hits fitting a track to all layers
0067         buildMatrixRecursive(m_cfg.n_layer-1,&hitMatrix,inputhits,outputTracks);
0068
0069     }
0070
0071
0072     void FarDetectorLinearTracking::buildMatrixRecursive(int level,
0073                                                         Eigen::MatrixXd* hitMatrix,
0074                                                         const std::vector<gsl::not_null<const edm4hep::TrackerHitCollection*>>& hits,
0075                                                         gsl::not_null<edm4eic::TrackSegmentCollection*> outputTracks ) const {
0076
0077       // Iterate over hits in this layer
0078       for(auto hit : (*hits[level])){
0079         auto pos = hit.getPosition();
0080         hitMatrix->col(level) << pos.x, pos.y, pos.z;
0081
0082         // Check the last two hits are within a certain angle of the optimum direction
0083         if(m_cfg.restrict_direction && level<m_cfg.n_layer-1){
0084           if(!checkHitPair(hitMatrix->col(level),hitMatrix->col(level+1))){
0085             continue;
0086           }
0087         }
0088
0089         if(level>0){
0090           buildMatrixRecursive(level-1,
0091                                hitMatrix,
0092                                hits,
0093                                outputTracks);
0094         }
0095         else{
0096           checkHitCombination(hitMatrix,outputTracks);
0097         }
0098       }
0099
0100     }
0101
0102
0103     void FarDetectorLinearTracking::checkHitCombination(Eigen::MatrixXd* hitMatrix,
0104                                                         gsl::not_null<edm4eic::TrackSegmentCollection*> outputTracks ) const {
0105
0106       Eigen::Vector3d weightedAnchor = (*hitMatrix)*m_layerWeights/(m_layerWeights.sum());
0107
0108       auto localMatrix = (*hitMatrix).colwise()-weightedAnchor;
0109
0110       Eigen::BDCSVD<Eigen::MatrixXd> svd(localMatrix.transpose(), Eigen::ComputeThinU | Eigen::ComputeThinV);
0111
0112       auto V = svd.matrixV();
0113
0114       // Rotate into principle components and calculate chi2/ndf
0115       auto rotatedMatrix = localMatrix.transpose()*V;
0116       auto residuals     = rotatedMatrix.rightCols(2);
0117       double chi2        = (residuals.array()*residuals.array()).sum()/(2*m_cfg.n_layer);
0118
0119       if(chi2>m_cfg.chi2_max) return;
0120
0121       edm4hep::Vector3d outPos = weightedAnchor.data();
0122       edm4hep::Vector3d outVec = V.col(0).data();
0123
0124       // Make sure fit was pointing in the right direction
0125       if(outVec.z>0) outVec = outVec*-1;
0126
0127       uint64_t          surface{0};         // Surface track was propagated to (possibly multiple per detector)
0128       uint32_t          system{0};          // Detector system track was propagated to
0129       edm4hep::Vector3f position(outPos.x,outPos.y,outPos.z);        // Position of the trajectory point [mm]
0130       edm4eic::Cov3f    positionError;      // Error on the position
0131       edm4hep::Vector3f momentum;       // 3-momentum at the point [GeV]
0132       edm4eic::Cov3f    momentumError;      // Error on the 3-momentum
0133       float             time{0};            // Time at this point [ns]
0134       float             timeError{0};       // Error on the time at this point
0135       float             theta = edm4eic::anglePolar(outVec);     // global polar direction of the fitted vector [rad]
0136       float             phi   = edm4eic::angleAzimuthal(outVec); // global azimuthal direction of the fitted vector [rad]
0137       edm4eic::Cov2f    directionError;     // Error on the polar and azimuthal angles
0138       float             pathlength{0};      // Pathlength from the origin to this point
0139       float             pathlengthError{0}; // Error on the pathlength
0140
0141       edm4eic::TrackPoint point({surface,system,position,positionError,momentum,momentumError,time,timeError,theta,phi,directionError,pathlength,pathlengthError});
0142
0143       float length      = 0;
0144       float lengthError = 0;
0145       auto segment = (*outputTracks)->create(length,lengthError);
0146
0147       segment.addToPoints(point);
0148
0149
0150     }
0151
0152     // Check if a pair of hits lies close to the optimum direction
0153     bool FarDetectorLinearTracking::checkHitPair(const Eigen::Vector3d& hit1,
0154                                                 const Eigen::Vector3d& hit2) const {
0155
0156       Eigen::Vector3d hitDiff = hit2-hit1;
0157       hitDiff.normalize();
0158
0159       double angle = std::acos(hitDiff.dot(m_optimumDirection));
0160
0161       debug("Vector: x={}, y={}, z={}",hitDiff.x(),hitDiff.y(),hitDiff.z());
0162       debug("Optimum: x={}, y={}, z={}",m_optimumDirection.x(),m_optimumDirection.y(),m_optimumDirection.z());
0163       debug("Angle: {}, Tolerance {}",angle,m_cfg.step_angle_tolerance);
0164
0165       if(angle>m_cfg.step_angle_tolerance) return false;
0166
0167       return true;
0168
0169     }
0170
0171 }