GSFProcessor.h

#pragma once
 
//--- Framework ---//
#include "Framework/Configure/Parameters.h"
#include "Framework/EventProcessor.h"
#include "Framework/Logger.h"
#include "Framework/RandomNumberSeedService.h"
 
//--- C++ ---//
#include <memory>
#include <random>
 
//--- LDMX ---//
#include "Tracking/Reco/TrackingGeometryUser.h"
 
//--- ACTS ---//
 
// Utils and Definitions
#include "Acts/Definitions/Common.hpp"
#include "Acts/Definitions/TrackParametrization.hpp"
#include "Acts/Definitions/Units.hpp"
#include "Acts/EventData/BoundTrackParameters.hpp"
#include "Acts/Utilities/Logger.hpp"
 
// geometry
#include "Acts/Geometry/GeometryContext.hpp"
 
// magfield
#include "Acts/MagneticField/MagneticFieldContext.hpp"
#include "Acts/MagneticField/MagneticFieldProvider.hpp"
 
// geometry
#include <Acts/Geometry/TrackingGeometry.hpp>
 
// propagation testing
#include "Acts/MagneticField/ConstantBField.hpp"
#include "Acts/Propagator/ActorList.hpp"
#include "Acts/Propagator/EigenStepperDenseExtension.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"
#include "Acts/Surfaces/PerigeeSurface.hpp"
#include "Acts/Utilities/Logger.hpp"
#include "Tracking/EigenStepper.h"  // IWYU pragma: keep
 
// Kalman Filter
 
// #include "Acts/EventData/Measurement.hpp"
#include "Acts/EventData/MultiTrajectory.hpp"
#include "Acts/EventData/MultiTrajectoryHelpers.hpp"
#include "Acts/EventData/VectorTrackContainer.hpp"
#include "Acts/Geometry/GeometryIdentifier.hpp"
#include "Acts/TrackFinding/CombinatorialKalmanFilter.hpp"
#include "Acts/TrackFinding/MeasurementSelector.hpp"
#include "Acts/TrackFitting/GainMatrixUpdater.hpp"
#include "Acts/Utilities/CalibrationContext.hpp"
 
//--- Refit with backward propagation ---//
#include "Acts/TrackFitting/KalmanFitter.hpp"
 
// GSF
#include "Acts/Propagator/MultiEigenStepperLoop.hpp"
#include "Acts/TrackFitting/BetheHeitlerApprox.hpp"
#include "Acts/TrackFitting/GaussianSumFitter.hpp"
#include "Acts/TrackFitting/GsfMixtureReduction.hpp"
 
//--- Tracking ---//
#include "Tracking/Event/Measurement.h"
#include "Tracking/Event/Track.h"
#include "Tracking/Reco/TrackExtrapolatorTool.h"
#include "Tracking/Sim/IndexSourceLink.h"
#include "Tracking/Sim/MeasurementCalibrator.h"
#include "Tracking/Sim/TrackingUtils.h"
 
//--- Interpolated magnetic field ---//
#include "Tracking/Sim/BFieldXYZUtils.h"
 
using ActionList =
    Acts::ActorList<Acts::detail::SteppingLogger, Acts::MaterialInteractor,
                    Acts::EndOfWorldReached>;
 
// using GsfPropagator = Acts::Propagator<
//                         Acts::MultiEigenStepperLoop<
//                           Acts::StepperExtensionList<
//                            Acts::detail::GenericDefaultExtension<double> >,
//                           Acts::WeightedComponentReducerLoop,
//                           Acts::detail::VoidAuctioneer>,
//                         Acts::Navigator>;
 
using MultiStepper = Acts::MultiEigenStepperLoop<>;
using Propagator = Acts::Propagator<Acts::EigenStepper<>, Acts::Navigator>;
using GsfPropagator = Acts::Propagator<MultiStepper, Acts::Navigator>;
using GsfExtrapPropagator =
    Acts::Propagator<Acts::EigenStepper<>, Acts::VoidNavigator>;
using GsfFitter =
    Acts::GaussianSumFitter<GsfPropagator, Acts::VectorMultiTrajectory>;
 
namespace tracking {
namespace reco {
 
class GSFProcessor final : public TrackingGeometryUser {
 public:
  GSFProcessor(const std::string& name, framework::Process& process);
 
  virtual ~GSFProcessor() = 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:
  // Forms the layer to acts map
  // auto makeLayerSurfacesMap(std::shared_ptr<const Acts::TrackingGeometry>
  // trackingGeometry) const -> std::unordered_map<unsigned int, const
  // Acts::Surface*>;
 
  // Make geoid -> source link map Measurements
  // auto makeGeoIdSourceLinkMap(
  //     const geo::TrackersTrackingGeometry& tg,
  //     const std::vector<ldmx::Measurement > &ldmxsps) ->
  //     std::unordered_multimap<Acts::GeometryIdentifier,
  //     ActsExamples::IndexSourceLink>;
 
  // If we want to dump the tracking geometry
  // bool dumpobj_{false};
  // int pionstates_{10};
  // int nevents_{0};
  // Processing time counter
  // double processing_time_{0.};
 
  int nevents_{0};
  int n_input_tracks_{0};
  int n_gsf_failed_{0};
  int n_output_tracks_{0};
  int n_target_extrap_failed_{0};
  int n_ecal_extrap_failed_{0};
  double processing_time_{0.};
 
  // Fallback statistics
  int n_fieldmap_gsf_failed_{0};
  int n_fallback_gsf_recovered_{0};
  int n_fieldmap_start_extrap_failed_{0};
  int n_fallback_start_extrap_recovered_{0};
  int n_fieldmap_target_extrap_failed_{0};
  int n_fallback_target_extrap_recovered_{0};
  int n_fieldmap_ecal_extrap_failed_{0};
  int n_fallback_ecal_extrap_recovered_{0};
  int n_start_extrap_failed_{0};
 
  // refitting of tracks
  // bool kf_refit_{false};
  // bool gsf_refit_{false};
 
  bool debug_{false};
 
  std::default_random_engine generator_;
 
  std::shared_ptr<std::normal_distribution<float>> normal_;
 
  // 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 use1Dmeasurements_{true};
 
  // Minimum number of hits on tracks
  // int min_hits_{7};
 
  std::vector<double> extrapolate_location_{0., 0., 0.};
  // bool use_seed_perigee_{false};
 
  std::string measurement_collection_{"TaggerMeasurements"};
 
  // double outlier_pval_{3.84};
 
  std::string out_trk_collection_{"GSFTracks"};
 
  // Select the hits using TrackID and pdg_id__
 
  // int track_id_{-1};
  // int pdg_id_{11};
 
  // Mass for the propagator hypothesis in MeV
  // double mass_{0.511};
 
  std::string seed_coll_name_{"seedTracks"};
 
  std::unique_ptr<const GsfFitter> gsf_;
 
  std::unique_ptr<const GsfFitter> gsf_zero_b_;
 
  std::unique_ptr<const GsfFitter> gsf_const_b_;
 
  std::string track_collection_{"TaggerTracks"};
 
  std::string meas_collection_{"DigiTaggerSimHits"};
 
  std::string track_passname_;
 
  std::string meas_passname_;
 
  std::string track_collection_event_passname_;
 
  std::string meas_collection_event_passname_;
 
  size_t max_components_{4};
 
  bool abort_on_error_{false};
 
  bool disable_all_material_handling_{false};
 
  double weight_cutoff_{1.0e-4};
 
  double propagator_step_size_{200.};
 
  int propagator_max_steps_{1000};
 
  std::string field_map_{""};
 
  double bfield_{-1.5};
  
  BFieldDistortion bfield_distortion_{};
 
  bool use_perigee_{false};
 
  // Keep track on which system this processor is running on
  bool tagger_tracking_{true};
 
  double tagger_start_x_{-617.};
 
  std::unique_ptr<const Propagator> propagator_;
 
  std::unordered_map<unsigned int, const Acts::Surface*> layer_surface_map_;
 
  // Track Extrapolator Tool (VoidNavigator to reach surfaces outside geometry)
  std::unique_ptr<const GsfExtrapPropagator> propagator_extrap_;
  std::shared_ptr<tracking::reco::TrackExtrapolatorTool<GsfExtrapPropagator>>
      trk_extrap_;
 
  // Zero-field extrapolator, fallback for the recoil
  std::unique_ptr<const GsfExtrapPropagator> propagator_extrap_zero_b_;
  std::shared_ptr<tracking::reco::TrackExtrapolatorTool<GsfExtrapPropagator>>
      trk_extrap_zero_b_;
 
  // Constant-field extrapolator, fallback for the tagger
  std::unique_ptr<const GsfExtrapPropagator> propagator_extrap_const_b_;
  std::shared_ptr<tracking::reco::TrackExtrapolatorTool<GsfExtrapPropagator>>
      trk_extrap_const_b_;
 
  std::shared_ptr<Acts::Surface> beam_origin_surface_;
 
  std::shared_ptr<Acts::Surface> tagger_start_surface_;
 
  std::shared_ptr<Acts::Surface> target_surface_;
 
  std::shared_ptr<Acts::Surface> ecal_surface_;
 
};  // GSFProcessor
 
}  // namespace reco
}  // namespace tracking