added AlignSegmentExtracter: this can be used to extract T segments from long...

added AlignSegmentExtracter: this can be used to extract T segments from long tracks. currently only works for PrKalmanFilter

Alignment/AlignTrTools/CMakeLists.txt

@@ -28,6 +28,7 @@ gaudi_add_module(AlignTrTools
 	src/AlignTrackCloneRemover.cpp
         src/ResettingSink.cpp
 	src/AlignVertexTrackExtracter.cpp
+	src/AlignSegmentExtracter.cpp
     LINK
         Alignment::AlignmentInterfacesLib
         Boost::headers

Alignment/AlignTrTools/src/AlignSegmentExtracter.cpp

+/*****************************************************************************\
+* (c) Copyright 2000-2018 CERN for the benefit of the LHCb Collaboration      *
+*                                                                             *
+* This software is distributed under the terms of the GNU General Public      *
+* Licence version 3 (GPL Version 3), copied verbatim in the file "COPYING".   *
+*                                                                             *
+* In applying this licence, CERN does not waive the privileges and immunities *
+* granted to it by virtue of its status as an Intergovernmental Organization  *
+* or submit itself to any jurisdiction.                                       *
+\*****************************************************************************/
+/** @class AlignTrackCloneRemover AlignTrackCloneRemover.h
+ *
+ *  Merge different track lists.
+ *
+ *  @author Wouter Hulsbergen
+ *  @date   16/03/2025
+ */
+
+#include "Event/Track.h"
+//#include "Event/StateParameters.h"
+#include "LHCbAlgs/Transformer.h"
+#include "PrKalmanFilter/KF.h"
+
+namespace LHCb {
+
+  struct AlignSegmentExtracter final : public Algorithm::Transformer<LHCb::Tracks( const LHCb::Track::Range& )> {
+
+    AlignSegmentExtracter( const std::string& name, ISvcLocator* pSvcLocator )
+        : Transformer( name, pSvcLocator, {"InputLocation", {}}, {"OutputLocation", {}} ) {}
+
+    LHCb::Tracks operator()( const LHCb::Track::Range& in ) const override {
+
+      LHCb::Tracks out;
+      out.reserve( in.size() );
+      for ( const auto& trkin : in ) {
+        // for now we assume that the input track is a fitted long track and that we want to extract the T segment, and
+        // refit only that part. let's also assume that it was fitted with PrKalmanFilter: then we can extract the nodes
+        // and to the work with that. that's just a lot more efficient than the stuff below.
+
+        // we can refine later.
+        auto prkfr = dynamic_cast<const LHCb::PrKalmanFitResult*>( trkin->fitResult() );
+        if ( prkfr ) {
+          // find the nodes corresponding to the first and last T hit
+          auto begin = prkfr->fitnodes.begin();
+          auto end   = prkfr->fitnodes.end();
+          --end;
+          while ( begin != prkfr->fitnodes.end() && !begin->lhcbID.isFT() ) ++begin;
+          while ( end != begin && !end->lhcbID.isFT() ) --end;
+          ++end;
+          // create a new fitresult
+          auto  fitresult = std::make_unique<LHCb::PrKalmanFitResult>( *prkfr );
+          auto& fitnodes  = fitresult->fitnodes;
+          fitnodes.clear();
+          fitnodes.insert( fitnodes.end(), begin, end );
+          fitresult->gain_matrices.clear();
+          fitresult->gain_matrices.resize( fitresult->fitnodes.size() );
+          fitresult->number_of_iter = -1;
+          // relabel all hits to be not outliers
+          for ( auto& node : fitnodes ) node.m_is_outlier = false;
+          // call the fit, but don't update transport. there is no function for hit in prkf so we copy what we need
+          LHCb::Pr::Tracks::Fit::KF::FitConfiguration fit_config{1.e2, 1.e2, 0.01, 0.01, 1e-4, 9999., 9., 9999., 1, 10};
+          const int                                   n_track_parameters = 5;
+          // This is silly: pre_fit_init does not set the off-diagonal elements. that's not very neat.
+          fitnodes.front().predicted_state_cov[LHCb::Pr::Tracks::Fit::Node::forward] = Gaudi::TrackSymMatrix{};
+          fitnodes.back().predicted_state_cov[LHCb::Pr::Tracks::Fit::Node::backward] = Gaudi::TrackSymMatrix{};
+          LHCb::Pr::Tracks::Fit::KF::pre_fit_init( fitnodes, fit_config );
+          auto       chi2      = LHCb::Pr::Tracks::Fit::KF::fit( fitnodes, n_track_parameters );
+          auto const tolerance = chi2.nDoF() * 0.01;
+          int        iter{2};
+          for ( bool has_converged{false}; iter <= fit_config.max_fit_iter && !has_converged; ++iter ) {
+            LHCb::Pr::Tracks::Fit::KF::smooth_and_update_ref_vector( fitnodes );
+            for ( auto& node : fitnodes ) node.project_reference();
+            auto const thischi2 = LHCb::Pr::Tracks::Fit::KF::fit( fitnodes, n_track_parameters );
+            auto const dchi2    = chi2.chi2() - thischi2.chi2();
+            chi2                = thischi2;
+            has_converged       = std::abs( dchi2 ) < tolerance;
+          }
+          fitresult->number_of_iter = iter;
+          // Remove outliers (above we specify at most one outlier)
+          int outlier_iter_buffer{0}, chi2_cut_buffer{0};
+          LHCb::Pr::Tracks::Fit::KF::remove_outliers( fitnodes, fit_config, chi2, outlier_iter_buffer, chi2_cut_buffer,
+                                                      n_track_parameters );
+          // Call the classical smoother
+          LHCb::Pr::Tracks::Fit::KF::classical_smoother_iteration_for_alignment( fitnodes, fitresult->gain_matrices );
+          // Create the new track
+          auto                      trkout = new LHCb::Track{trkin->history(), LHCb::Event::Enum::Track::Type::Ttrack,
+                                        LHCb::Event::Enum::Track::PatRecStatus::PatRecIDs};
+          std::vector<LHCb::LHCbID> lhcbids;
+          // we could copy these from the nodes, but from the track is more efficient because already ordered
+          for ( const auto lhcbid : trkin->lhcbIDs() )
+            if ( lhcbid.isFT() ) lhcbids.emplace_back( lhcbid );
+          // create a new track
+          trkout->setLhcbIDs( lhcbids );
+          auto stateA = state( fitnodes.front() );
+          auto stateB = state( fitnodes.back() );
+          if ( stateA.z() > stateB.z() ) std::swap( stateA, stateB );
+          trkout->addToStates( stateA );
+          trkout->addToStates( stateB );
+          trkout->setChi2AndDoF( chi2.chi2(), chi2.nDoF() );
+          trkout->setFitResult( fitresult.release() );
+          trkout->setFitStatus( LHCb::Event::Enum::Track::FitStatus::Fitted );
+          out.insert( trkout );
+        }
+      }
+      numtracksin += in.size();
+      numtracksout += out.size();
+      return out;
+    }
+
+    // Keep some counters
+    mutable Gaudi::Accumulators::StatCounter<> numtracksin{this, "Number of input tracks"};
+    mutable Gaudi::Accumulators::StatCounter<> numtracksout{this, "Number of output segments"};
+  };
+
+  DECLARE_COMPONENT_WITH_ID( AlignSegmentExtracter, "AlignTTrackSegmentExtracter" )
+} // namespace LHCb