TrkrClusterOutlierFittingAlgorithm.hpp

Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file TrkrClusterOutlierFittingAlgorithm.hpp

#pragma once

#include <iostream>
#include <map>
#include <random>
#include <stdexcept>
#include <boost/program_options.hpp>

#include "Acts/TrackFitting/GainMatrixSmoother.hpp"
#include "Acts/TrackFitting/GainMatrixUpdater.hpp"
#include "Acts/Geometry/GeometryIdentifier.hpp"
#include "Acts/MagneticField/ConstantBField.hpp"
#include "Acts/MagneticField/InterpolatedBFieldMap.hpp"
#include "Acts/MagneticField/SharedBField.hpp"
#include "Acts/Propagator/EigenStepper.hpp"
#include "Acts/Propagator/Navigator.hpp"
#include "Acts/Propagator/Propagator.hpp"
#include "Acts/Surfaces/Surface.hpp"
#include "Acts/Utilities/Helpers.hpp"
#include "Acts/Utilities/ParameterDefinitions.hpp"
#include "Acts/TrackFitting/KalmanFitter.hpp"
#include "Acts/Geometry/TrackingGeometry.hpp"

#include "ActsExamples/Plugins/BField/ScalableBField.hpp"
#include "ActsExamples/EventData/Track.hpp"
#include "ActsExamples/Framework/BareAlgorithm.hpp"
#include "ActsExamples/Plugins/BField/BFieldOptions.hpp"
#include "ActsExamples/EventData/TrkrClusterSourceLink.hpp"

struct ResidualOutlierFinder {
  std::vector<std::pair<long unsigned int, double>> chi2Cuts;
  int verbosity = 0;

  template <typename track_state_t>
  bool operator()(const track_state_t& state) const {
    
    if(not state.hasCalibrated() or not state.hasPredicted())
      return false;

    // residuals = state.calibrated() - state.projector() * state.predicted();
    auto distance = (state.calibrated() - state.projector() * state.predicted()).norm();

    const auto& projector = state.projector();
    const auto& predictedCov = state.predictedCovariance();
    const auto& stateCov = state.effectiveCalibratedCovariance();
    
    double chi2 = ((state.calibrated() - state.projector() * state.predicted()).transpose() * 
       (( stateCov + projector * predictedCov * projector.transpose())).inverse() * (state.calibrated() - state.projector() * state.predicted()))
                   .eval()(0,0);

    if(verbosity > 2)
      {
  std::cout << "geoid : " << state.referenceSurface().geometryId() << std::endl;
  std::cout << "Distance between prediction and measurement is: " << distance << std::endl;
  std::cout << "chi2 " << chi2  << std::endl;
      }

    auto volID = state.referenceSurface().geometryId().volume();

    for(auto& pair : chi2Cuts)
      {
  if(verbosity > 2)
    {
      std::cout << "VolID : " << volID << " checking against " << pair.first 
          << ",   chicut " << pair.second << " with chi2 " << chi2 << std::endl;
    }

  if(volID == pair.first)
    {
      return chi2 >= pair.second;
    }
      }
    
    return false;
  }
};


namespace ActsExamples {

class TrkrClusterOutlierFittingAlgorithm : public BareAlgorithm
{
 public:
  using FitterResult
    = Acts::Result<Acts::KalmanFitterResult<ActsExamples::TrkrClusterSourceLink>>;
  using FitterFunction
    = std::function<FitterResult(
        const std::vector<ActsExamples::TrkrClusterSourceLink>&,
        const TrackParameters&,
        const Acts::KalmanFitterOptions<ResidualOutlierFinder>&)>;

  using DirectedFitterFunction
    = std::function<FitterResult(
        const std::vector<ActsExamples::TrkrClusterSourceLink>&,
        const TrackParameters&,
        const Acts::KalmanFitterOptions<ResidualOutlierFinder>&,
        const std::vector<const Acts::Surface*>&)>;
      

  static FitterFunction makeFitterFunction(
      std::shared_ptr<const Acts::TrackingGeometry> trackingGeometry,
      Options::BFieldVariant magneticField);

  static DirectedFitterFunction makeFitterFunction(
      Options::BFieldVariant magneticField);

  struct Config
  {
    FitterFunction fit;
    DirectedFitterFunction dFit;
  };

  TrkrClusterOutlierFittingAlgorithm(Config cfg, Acts::Logging::Level lvl);


private:
  Config m_cfg;
};

}