CKFProcessor.h

#pragma once
 
//--- Framework ---//
#include "Framework/Configure/Parameters.h"
#include "Framework/EventProcessor.h"
#include "Framework/Logger.h"
 
//--- C++ ---//
#include <map>
#include <memory>
 
//--- LDMX ---//
#include "Tracking/EigenStepper.h"  // IWYU pragma: keep
#include "Tracking/Reco/TrackingGeometryUser.h"
 
//--- ACTS ---//
 
// Utils and Definitions
 
// geometry
 
// magfield
 
// geometry
 
// propagation testing
#include "Acts/Propagator/ActorList.hpp"
#include "Acts/Propagator/MaterialInteractor.hpp"
#include "Acts/Propagator/Navigator.hpp"
#include "Acts/Propagator/Propagator.hpp"
#include "Acts/Propagator/StandardAborters.hpp"
#include "Acts/Propagator/VoidNavigator.hpp"
#include "Acts/Propagator/detail/SteppingLogger.hpp"
 
// Kalman Filter
 
// #include "Acts/EventData/Measurement.hpp"
#include "Acts/EventData/VectorTrackContainer.hpp"
#include "Acts/Geometry/GeometryIdentifier.hpp"
#include "Acts/TrackFinding/CombinatorialKalmanFilter.hpp"
#include "Acts/TrackFitting/GainMatrixUpdater.hpp"
 
//--- Refit with backward propagation ---//
 
// GSF
// #include "Acts/TrackFitting/GaussianSumFitter.hpp"
 
//--- Tracking ---//
#include "Tracking/Event/Measurement.h"
#include "Tracking/Event/Track.h"
#include "Tracking/Reco/TrackExtrapolatorTool.h"
#include "Tracking/Sim/IndexSourceLink.h"
 
//--- Interpolated magnetic field ---//
#include "Tracking/Sim/BFieldXYZUtils.h"
// mg Aug 2024 not sure if these are needed...
using Updater = Acts::GainMatrixUpdater;
 
using ActionList =
    Acts::ActorList<Acts::detail::SteppingLogger, Acts::MaterialInteractor,
                    Acts::EndOfWorldReached>;
 
using CkfPropagator = Acts::Propagator<Acts::EigenStepper<>, Acts::Navigator>;
using ExtrapPropagator =
    Acts::Propagator<Acts::EigenStepper<>, Acts::VoidNavigator>;
using TrackContainer = Acts::TrackContainer<Acts::VectorTrackContainer,
                                            Acts::VectorMultiTrajectory>;
 
namespace tracking {
namespace reco {
 
class CKFProcessor final : public TrackingGeometryUser {
 public:
  CKFProcessor(const std::string& name, framework::Process& process);
 
  virtual ~CKFProcessor() = default;
 
  void onProcessStart() override;
 
  void onNewRun(const ldmx::RunHeader& rh) override;
 
  void onProcessEnd() override;
 
  void configure(framework::config::Parameters& parameters) override;
 
  void produce(framework::Event& event) override;
 
 private:
  // Make geoid -> source link map Measurements
  auto makeGeoIdSourceLinkMap(const geo::TrackersTrackingGeometry& tg,
                              const std::vector<ldmx::Measurement>& ldmxsps)
      -> std::unordered_multimap<Acts::GeometryIdentifier,
                                 acts_examples::IndexSourceLink>;
 
  template <typename geometry_t, typename source_link_hash_t,
            typename source_link_equality_t>
  std::vector<std::vector<std::size_t>> computeSharedHits(
      std::vector<ldmx::Track> tracks, std::vector<ldmx::Measurement> meas_coll,
      geometry_t& tg, source_link_hash_t&& sourceLinkHash,
      source_link_equality_t&& sourceLinkEquality) const;
 
  // If we want to dump the tracking geometry
  bool dumpobj_{false};
 
  int pionstates_{10};
 
  int nevents_{0};
 
  // Processing time counter
  double processing_time_{0.};
 
  // time profiling
  std::map<std::string, double> profiling_map_;
 
  bool debug_acts_{false};
 
  std::shared_ptr<Acts::PlaneSurface> target_surface_;
  Acts::RotationMatrix3 surf_rotation_;
  // Constant BField
  double bfield_{0};
  // Use constant bfield
  bool const_b_field_{true};
 
  // Remove stereo measurements
  bool remove_stereo_{false};
 
  // Use 2d measurements instead of 1D
  bool use1_dmeasurements_{true};
 
  // Minimum number of hits on tracks
  int min_hits_{7};
 
  // Stepping size (in mm)
  double propagator_step_size_{200.};
  int propagator_max_steps_{1000};
 
  // The extrapolation surface
  bool use_extrapolate_location_{true};
  std::vector<double> extrapolate_location_{0., 0., 0.};
  bool use_seed_perigee_{false};
 
  // The measurement collection to use for track reconstruction
  std::string measurement_collection_{"TaggerMeasurements"};
 
  std::string sim_particles_coll_name_;
  std::string sim_particles_event_passname_;
 
  // Outlier removal pvalue
  // The Chi2Cut is applied at filtering stage.
  // 1DOF pvalues: 0.1 = 2.706 0.05 = 3.841 0.025 = 5.024 0.01 = 6.635 0.005
  // = 7.879 The probability to reject a good measurement is pvalue The
  // probability to reject an outlier is given in NIM A262 (1987) 444-450
  double outlier_pval_;
 
  // The output track collection
  std::string out_trk_collection_{"Tracks"};
 
  // The seed track collection
  std::string seed_coll_name_{"seedTracks"};
 
  // The interpolated bfield
  std::string field_map_{""};
 
  std::string input_pass_name_{""};
 
  // The Propagator
  std::unique_ptr<const CkfPropagator> propagator_;
 
  // The CKF
  std::unique_ptr<
      const Acts::CombinatorialKalmanFilter<CkfPropagator, TrackContainer>>
      ckf_;
 
  // Track Extrapolator Tool (uses VoidNavigator to propagate freely to any
  // surface)
  std::unique_ptr<const ExtrapPropagator> propagator_extrap_;
  std::shared_ptr<tracking::reco::TrackExtrapolatorTool<ExtrapPropagator>>
      trk_extrap_;
 
  // Zero-B CKF as fallback
  std::unique_ptr<const CkfPropagator> propagator_zero_b_;
  std::unique_ptr<
      const Acts::CombinatorialKalmanFilter<CkfPropagator, TrackContainer>>
      ckf_zero_b_;
  std::unique_ptr<const ExtrapPropagator> propagator_extrap_zero_b_;
  std::shared_ptr<tracking::reco::TrackExtrapolatorTool<ExtrapPropagator>>
      trk_extrap_zero_b_;
 
  // Const-B (1.5T) CKF as fallback for tagger
  std::unique_ptr<const CkfPropagator> propagator_const_b_;
  std::unique_ptr<
      const Acts::CombinatorialKalmanFilter<CkfPropagator, TrackContainer>>
      ckf_const_b_;
  std::unique_ptr<const ExtrapPropagator> propagator_extrap_const_b_;
  std::shared_ptr<tracking::reco::TrackExtrapolatorTool<ExtrapPropagator>>
      trk_extrap_const_b_;
 
  int nseeds_;
  int ntracks_;
  int eventnr_;
 
  // CKF fallback statistics
  int n_fieldmap_ckf_failed_tagger_{0};
  int n_fieldmap_ckf_failed_recoil_{0};
  int n_constb_ckf_recovered_tagger_{0};
  int n_zerob_ckf_recovered_recoil_{0};
 
  // Extrapolation fallback statistics
  int n_fieldmap_target_extrap_failed_tagger_{0};
  int n_constb_target_extrap_recovered_tagger_{0};
  int n_fieldmap_target_extrap_failed_recoil_{0};
  int n_zerob_target_extrap_recovered_recoil_{0};
  int n_fieldmap_ecal_extrap_failed_recoil_{0};
  int n_zerob_ecal_extrap_recovered_recoil_{0};
 
  // BField Systematics
  BFieldDistortion bfield_distortion_{};
 
  // Keep track on which system this processor is running on
  bool tagger_tracking_{true};
 
  // Target point constraint (recoil only), same logic as the seeder
  bool use_target_constraint_{false};
  // extra target point sigma [mm] and its 1/p part [mm GeV]
  std::vector<double> target_sigma_extra_{0., 0.};
  std::vector<double> target_sigma_extra_ms_{0., 0.};
  // target thickness for the vertex depth term [mm]
  double target_thickness_{0.};
  // apply the target point only up to this many hits, -1 for all
  int target_constraint_max_hits_{-1};
  std::string tagger_trks_collection_{""};
  std::string tagger_trks_event_collection_passname_{""};
 
  struct TargetPoint {
    Acts::Vector2 loc_;
    Acts::Matrix<2, 2> cov_;
  };
 
  std::vector<TargetPoint> targetPoints(framework::Event& event);
 
  // Target constraint statistics
  int n_tgt_tagger_{0};
  int n_tgt_nocov_{0};
 
  // zero field outside the map instead of an out of bounds error
  bool field_zero_outside_{false};
 
  // field-map CKF failures by error message
  std::map<std::string, int> ckf_errors_;
 
};  // CKFProcessor
 
}  // namespace reco
}  // namespace tracking