GlobalChi2Fitter.h

Go to the documentation of this file.

/*
  Copyright (C) 2002-2024 CERN for the benefit of the ATLAS collaboration
*/
 
#ifndef GLOBALCHI2FITTER_H
#define GLOBALCHI2FITTER_H
//#define GXFDEBUGCODE
#include "TrkDetDescrInterfaces/IMaterialEffectsOnTrackProvider.h"
#include "AthenaBaseComps/AthAlgTool.h"
#include "AthenaBaseComps/AthCheckedComponent.h"
#include "GaudiKernel/ToolHandle.h"
#include "GaudiKernel/EventContext.h"
 
 
#include "TrkToolInterfaces/ITrkMaterialProviderTool.h"
#include "TrkToolInterfaces/IResidualPullCalculator.h"
#include "TrkToolInterfaces/IRIO_OnTrackCreator.h"
#include "TrkToolInterfaces/IUpdator.h"
#include "TrkToolInterfaces/IBoundaryCheckTool.h"
 
#include "TrkExInterfaces/IExtrapolator.h"
#include "TrkExInterfaces/IPropagator.h"
#include "TrkExInterfaces/INavigator.h"
#include "TrkExInterfaces/IMultipleScatteringUpdator.h"
#include "TrkExInterfaces/IEnergyLossUpdator.h"
#include "TrkExInterfaces/IMaterialEffectsUpdator.h"
 
#include "TrkFitterInterfaces/IGlobalTrackFitter.h"
 
#include "TrkGlobalChi2Fitter/GXFTrajectory.h"
#include "TrkMaterialOnTrack/MaterialEffectsOnTrack.h"
#include "TrkFitterUtils/FitterStatusCode.h"
#include "TrkEventPrimitives/PropDirection.h"
#include "MagFieldConditions/AtlasFieldCacheCondObj.h"
#include "MagFieldElements/AtlasFieldCache.h"
 
#include "StoreGate/ReadCondHandleKey.h"
#include "TrkGeometry/TrackingGeometry.h"
 
#include "TrkEventUtils/ClusterSplitProbabilityContainer.h"
#include "TrkEventPrimitives/TransportJacobian.h"
 
#include <memory>
#include <mutex>
 
class AtlasDetectorID;
 
 
namespace Trk {
  class Track;
  class IMagneticFieldTool;
  class MeasuredPerigee;
  class PrepRawDataComparisonFunction;
  class MeasurementBaseComparisonFunction;
  class MaterialEffectsOnTrack;
  class Layer;
  class CylinderLayer;
  class DiscLayer;
  class MagneticFieldProperties;
  class TrackingVolume;
  class Volume;
 
  class GlobalChi2Fitter: public extends<AthCheckedComponent<AthAlgTool>, IGlobalTrackFitter> {
    struct PropagationResult {
      std::unique_ptr<const TrackParameters> m_parameters;
      std::optional<TransportJacobian> m_jacobian;
      std::optional<std::vector<std::unique_ptr<TrackParameters>>> m_preholes;
    };
 
    /*
     * This struct serves as a very simple container for the five different
     * hole and dead module states that we want to count, which are the dead
     * Pixels, dead SCTs, Pixel holes, SCT holes, and SCT double holes.
     */
    struct TrackHoleCount {
      unsigned int m_pixel_hole = 0;
      unsigned int m_sct_hole = 0;
      unsigned int m_sct_double_hole = 0;
      unsigned int m_pixel_dead = 0;
      unsigned int m_sct_dead = 0;
    };
 
    enum FitterStatusType {
      S_FITS,
      S_SUCCESSFUL_FITS,
      S_MAT_INV_FAIL,
      S_NOT_ENOUGH_MEAS,
      S_PROPAGATION_FAIL,
      S_INVALID_ANGLES,
      S_NOT_CONVERGENT,
      S_HIGH_CHI2,
      S_LOW_MOMENTUM,
      __S_MAX_VALUE
    };
 
    struct Cache {
      /*
       * Currently the information about what type of fit is being passed by the
       * presence of a TrackingVolume.
       */
      template <class T>
      static
      void objVectorDeleter(const std::vector<const T *> *ptr) {
        if (ptr) {
          for (const T *elm : *ptr) { delete elm; }
          delete ptr;
        }
      }
 
      const TrackingGeometry *m_trackingGeometry = nullptr;
      const TrackingVolume *m_caloEntrance = nullptr;
      const TrackingVolume *m_msEntrance = nullptr;
 
      bool m_calomat, m_extmat;
      bool m_idmat = true;
      bool m_sirecal;
      bool m_getmaterialfromtrack;
      bool m_reintoutl;
      bool m_matfilled = false;
      bool m_acceleration;
      bool m_fiteloss;
      bool m_asymeloss;
 
      std::vector<double> m_phiweight;
      std::vector<int> m_firstmeasurement;
      std::vector<int> m_lastmeasurement;
 
      std::vector < const Trk::Layer * >m_negdiscs;
      std::vector < const Trk::Layer * >m_posdiscs;
      std::vector < const Trk::Layer * >m_barrelcylinders;
 
      bool m_fastmat = true;
 
      int m_lastiter;
      int m_miniter;
 
      #ifdef GXFDEBUGCODE
      int m_iterations = 0;
      #endif
 
      Amg::MatrixX m_derivmat;
      Amg::SymMatrixX m_fullcovmat;
 
      std::vector< std::unique_ptr< const std::vector < const TrackStateOnSurface *>,
                                    void (*)(const std::vector<const TrackStateOnSurface *> *) > >
        m_matTempStore;
 
      MagField::AtlasFieldCache m_field_cache;
 
      FitterStatusCode m_fittercode;
 
      std::array<unsigned int, __S_MAX_VALUE> m_fit_status {};
      std::array<std::atomic<unsigned int>, __S_MAX_VALUE>  *m_fit_status_out = nullptr;
 
       Cache(const GlobalChi2Fitter *fitter):
        m_calomat(fitter->m_calomat),
        m_extmat(fitter->m_extmat),
        m_sirecal(fitter->m_sirecal),
        m_getmaterialfromtrack(fitter->m_getmaterialfromtrack),
        m_reintoutl(fitter->m_reintoutl),
        m_acceleration(fitter->m_acceleration),
        m_fiteloss(fitter->m_fiteloss),
        m_asymeloss(fitter->m_asymeloss),
        m_miniter(fitter->m_miniter),
        m_fit_status_out(&fitter->m_fit_status)
      {}
      ~Cache() {
         unsigned int idx=0;
         for (unsigned int a_fit_status : m_fit_status) {
            if (a_fit_status) {
               (*m_fit_status_out)[idx] += a_fit_status;
            }
            ++idx;
         }
      }
 
      Cache & operator=(const Cache &) = delete;
 
      void incrementFitStatus(enum FitterStatusType status)  {
         m_fit_status[status]++;
      }
    };
 
  public:
    GlobalChi2Fitter(
      const std::string &,
      const std::string &,
      const IInterface *
    );
 
    virtual ~ GlobalChi2Fitter();
 
    virtual StatusCode initialize() override;
    virtual StatusCode finalize() override;
 
    virtual std::unique_ptr<Track> fit(
      const EventContext& ctx,
      const PrepRawDataSet&,
      const TrackParameters&,
      const RunOutlierRemoval runOutlier = false,
      const ParticleHypothesis matEffects = nonInteracting
      ) const override final;
 
    virtual std::unique_ptr<Track> fit(
      const EventContext& ctx,
      const Track &,
      const RunOutlierRemoval runOutlier = false,
      const ParticleHypothesis matEffects = nonInteracting
    ) const override final;
 
    virtual std::unique_ptr<Track> fit(
      const EventContext& ctx,
      const MeasurementSet &,
      const TrackParameters &,
      const RunOutlierRemoval runOutlier = false,
      const ParticleHypothesis matEffects = nonInteracting
    ) const override final;
 
    virtual std::unique_ptr<Track> fit(
      const EventContext& ctx,
      const Track &,
      const PrepRawDataSet &,
      const RunOutlierRemoval runOutlier = false,
      const ParticleHypothesis matEffects = nonInteracting
    ) const override final;
 
    virtual std::unique_ptr<Track> fit(
      const EventContext& ctx,
      const Track &,
      const Track &,
      const RunOutlierRemoval runOutlier = false,
      const ParticleHypothesis matEffects = nonInteracting
    ) const override final;
 
    virtual std::unique_ptr<Track> fit(
      const EventContext& ctx,
      const Track &,
      const MeasurementSet &,
      const RunOutlierRemoval runOutlier = false,
      const ParticleHypothesis matEffects = nonInteracting
    ) const override final;
 
    virtual Track* alignmentFit(
      AlignmentCache&,
      const Track&,
      const RunOutlierRemoval  runOutlier=false,
      const ParticleHypothesis matEffects=Trk::nonInteracting
    ) const override;
 
  private:
 
    Track * fitIm(
      const EventContext& ctx,
      Cache & cache,
      const Track & inputTrack,
      const RunOutlierRemoval runOutlier,
      const ParticleHypothesis matEffects
    ) const;
 
    Track *myfit(
      const EventContext& ctx,
      Cache &,
      GXFTrajectory &,
      const TrackParameters &,
      const RunOutlierRemoval runOutlier = false,
      const ParticleHypothesis matEffects = nonInteracting
    ) const;
 
    Track *myfit_helper(
      Cache &,
      GXFTrajectory &,
      const TrackParameters &,
      const RunOutlierRemoval runOutlier = false,
      const ParticleHypothesis matEffects = nonInteracting
    ) const;
 
    Track *mainCombinationStrategy(
      const EventContext& ctx,
      Cache &,
      const Track &,
      const Track &,
      GXFTrajectory &,
      std::vector<MaterialEffectsOnTrack> &
    ) const;
 
    Track *backupCombinationStrategy(
      const EventContext& ctx,
      Cache &,
      const Track &,
      const Track &,
      GXFTrajectory &,
      std::vector<MaterialEffectsOnTrack> &
    ) const;
 
    void makeProtoState(
      Cache &,
      GXFTrajectory &,
      const TrackStateOnSurface *,
      int index = -1
     ) const;
 
    void makeProtoStateFromMeasurement(
      Cache &,
      GXFTrajectory &,
      const MeasurementBase *,
      const TrackParameters * trackpar = nullptr,
      bool isoutlier = false,
      int index = -1
    ) const;
 
    bool processTrkVolume(
      Cache &,
      const Trk::TrackingVolume * tvol
    ) const;
 
    static std::optional<std::pair<Amg::Vector3D, double>> addMaterialFindIntersectionDisc(
      Cache & cache,
      const DiscSurface & surface,
      const TrackParameters & param1,
      const TrackParameters & param2,
      const ParticleHypothesis mat
    ) ;
 
    static std::optional<std::pair<Amg::Vector3D, double>> addMaterialFindIntersectionCyl(
      Cache & cache,
      const CylinderSurface & surface,
      const TrackParameters & param1,
      const TrackParameters & param2,
      const ParticleHypothesis mat
    ) ;
 
    void addMaterialUpdateTrajectory(
      Cache & cache,
      GXFTrajectory & track,
      int offset,
      std::vector<std::pair<const Layer *, const Layer *>> & layers,
      const TrackParameters * ref1,
      const TrackParameters * ref2,
      ParticleHypothesis mat
    ) const;
 
    static void addMaterialGetLayers(
      Cache & cache,
      std::vector<std::pair<const Layer *, const Layer *>> & layers,
      std::vector<std::pair<const Layer *, const Layer *>> & uplayers,
      const std::vector<std::unique_ptr<GXFTrackState>> & states,
      GXFTrackState & first,
      GXFTrackState & last,
      const TrackParameters * refpar,
      bool hasmat
    ) ;
 
    void addIDMaterialFast(
      const EventContext& ctx,
      Cache & cache,
      GXFTrajectory & track,
      const TrackParameters * parameters,
      ParticleHypothesis part
    ) const;
 
    void addMaterial(
      const EventContext& ctx,
      Cache &,
      GXFTrajectory &,
      const TrackParameters *,
      ParticleHypothesis
    ) const;
 
    std::vector<std::unique_ptr<TrackParameters>> holesearchExtrapolation(
      const EventContext & ctx,
      const TrackParameters & src,
      const GXFTrackState & dst,
      PropDirection propdir
    ) const;
 
    std::unique_ptr<const TrackParameters> makePerigee(
      Cache &,
      const TrackParameters &,
      const ParticleHypothesis
    ) const;
 
    static void makeTrackFillDerivativeMatrix(
      Cache &,
      GXFTrajectory &
    ) ;
 
    std::unique_ptr<const TrackParameters> makeTrackFindPerigeeParameters(
      const EventContext &,
      Cache &,
      GXFTrajectory &,
      const ParticleHypothesis
    ) const;
 
    std::unique_ptr<GXFTrackState> makeTrackFindPerigee(
      const EventContext &,
      Cache &,
      GXFTrajectory &,
      const ParticleHypothesis
    ) const;
 
    std::unique_ptr<Track> makeTrack(
      const EventContext& ctx,
      Cache &,
      GXFTrajectory &,
      const ParticleHypothesis
    ) const;
 
    void fillResiduals(
      const EventContext& ctx,
      Cache &,
      GXFTrajectory &,
      int,
      Amg::SymMatrixX &,
      Amg::VectorX &,
      Amg::SymMatrixX &,
      bool &
    ) const;
 
    void fillDerivatives(
      GXFTrajectory & traj,
      bool onlybrem = false
    ) const;
 
    FitterStatusCode runIteration(
      const EventContext& ctx,
      Cache &,
      GXFTrajectory &,
      int,
      Amg::SymMatrixX &,
      Amg::VectorX &,
      Amg::SymMatrixX &,
      bool &
    ) const;
 
    FitterStatusCode updateFitParameters(
      GXFTrajectory &,
      Amg::VectorX &,
      const Amg::SymMatrixX &
    ) const;
 
 
    void updatePixelROTs(
      GXFTrajectory &,
      Amg::SymMatrixX &,
      Amg::VectorX &,
      const EventContext& evtctx
    ) const;
 
 
    GXFTrajectory *runTrackCleanerSilicon(
      const EventContext& ctx,
      Cache &,
      GXFTrajectory &,
      Amg::SymMatrixX &,
      Amg::SymMatrixX &,
      Amg::VectorX &,
      bool
    ) const;
 
    void runTrackCleanerTRT(
      Cache &,
      GXFTrajectory &,
      Amg::SymMatrixX &,
      Amg::VectorX &,
      Amg::SymMatrixX &,
      bool, bool, int,
      const EventContext& ctx
    ) const;
 
    PropagationResult calculateTrackParametersPropagateHelper(
      const EventContext &,
      const TrackParameters &,
      const GXFTrackState &,
      PropDirection,
      const MagneticFieldProperties&,
      bool,
      bool
    ) const;
 
    PropagationResult calculateTrackParametersPropagate(
      const EventContext &,
      const TrackParameters &,
      const GXFTrackState &,
      PropDirection,
      const MagneticFieldProperties&,
      bool,
      bool
    ) const;
 
    std::vector<std::reference_wrapper<GXFTrackState>> holeSearchStates(
      GXFTrajectory & trajectory
    ) const;
 
    std::optional<GlobalChi2Fitter::TrackHoleCount> holeSearchProcess(
      const EventContext & ctx,
      const std::vector<std::reference_wrapper<GXFTrackState>> & states
    ) const;
 
    void holeSearchHelper(
      const std::vector<std::unique_ptr<TrackParameters>> & hc,
      std::set<Identifier> & id_set,
      std::set<Identifier> & sct_set,
      TrackHoleCount & rv,
      bool count_holes,
      bool count_dead
    ) const;
 
    FitterStatusCode calculateTrackParameters(
      const EventContext& ctx,
      GXFTrajectory&,
      bool) const;
 
    std::variant<std::unique_ptr<const TrackParameters>, FitterStatusCode> updateEnergyLoss(
      const Surface &,
      const GXFMaterialEffects &,
      const TrackParameters &,
      double,
      int
    ) const;
 
    static void calculateDerivatives(GXFTrajectory &) ;
 
    void calculateTrackErrors(GXFTrajectory &, Amg::SymMatrixX &, bool) const;
 
    std::optional<TransportJacobian> numericalDerivatives(
      const EventContext& ctx,
      const TrackParameters *,
      const Surface &,
      PropDirection,
      const MagneticFieldProperties&
    ) const;
 
    virtual int iterationsOfLastFit() const;
 
    virtual void setMinIterations(int);
 
    static bool correctAngles(double &, double &) ;
 
    bool isMuonTrack(const Track &) const;
 
    void initFieldCache(
      const EventContext& ctx,
      Cache & cache
    ) const;
 
    ToolHandle<IRIO_OnTrackCreator> m_ROTcreator {this, "RotCreatorTool", "", ""};
    ToolHandle<IRIO_OnTrackCreator> m_broadROTcreator {this, "BroadRotCreatorTool", "", ""};
    ToolHandle<IUpdator> m_updator {this, "MeasurementUpdateTool", "", ""};
    ToolHandle<IExtrapolator> m_extrapolator {this, "ExtrapolationTool", "Trk::Extrapolator/CosmicsExtrapolator", ""};
    ToolHandle<IMultipleScatteringUpdator> m_scattool {this, "MultipleScatteringTool", "Trk::MultipleScatteringUpdator/AtlasMultipleScatteringUpdator", ""};
    ToolHandle<IEnergyLossUpdator> m_elosstool {this, "EnergyLossTool", "Trk::EnergyLossUpdator/AtlasEnergyLossUpdator", ""};
    ToolHandle<IMaterialEffectsUpdator> m_matupdator {this, "MaterialUpdateTool", "", ""};
    ToolHandle<IPropagator> m_propagator {this, "PropagatorTool", "", ""};
    ToolHandle<INavigator> m_navigator {this, "NavigatorTool", "Trk::Navigator/CosmicsNavigator", ""};
    ToolHandle<IResidualPullCalculator> m_residualPullCalculator {this, "ResidualPullCalculatorTool", "Trk::ResidualPullCalculator/ResidualPullCalculator", ""};
    ToolHandle<Trk::ITrkMaterialProviderTool> m_caloMaterialProvider {this, "CaloMaterialProvider", "Trk::TrkMaterialProviderTool/TrkMaterialProviderTool", ""};
    ToolHandle<IMaterialEffectsOnTrackProvider> m_calotool {this, "MuidTool", "Rec::MuidMaterialEffectsOnTrackProvider/MuidMaterialEffectsOnTrackProvider", ""};
    ToolHandle<IMaterialEffectsOnTrackProvider> m_calotoolparam {this, "MuidToolParam", "", ""};
    ToolHandle<IBoundaryCheckTool> m_boundaryCheckTool {this, "BoundaryCheckTool", "", "Boundary checking tool for detector sensitivities" };
 
    void throwFailedToGetTrackingGeomtry() const;
    const TrackingGeometry* trackingGeometry(Cache& cache,
                                             const EventContext& ctx) const
    {
      if (!cache.m_trackingGeometry)
        cache.m_trackingGeometry = retrieveTrackingGeometry(ctx);
      return cache.m_trackingGeometry;
    }
    const TrackingGeometry* retrieveTrackingGeometry(
      const EventContext& ctx) const
    {
      SG::ReadCondHandle<TrackingGeometry> handle(m_trackingGeometryReadKey,
                                                  ctx);
      if (!handle.isValid()) {
        throwFailedToGetTrackingGeomtry();
      }
      return handle.cptr();
    }
 
    SG::ReadCondHandleKey<TrackingGeometry> m_trackingGeometryReadKey{
      this,
      "TrackingGeometryReadKey",
      "AtlasTrackingGeometry",
      "Key of the TrackingGeometry conditions data."
    };
 
    SG::ReadCondHandleKey<AtlasFieldCacheCondObj> m_field_cache_key{
      this,
      "AtlasFieldCacheCondObj",
      "fieldCondObj",
      "Trk::GlobalChi2Fitter field conditions object key"
    };
 
    const AtlasDetectorID *m_DetID = nullptr;
 
    Gaudi::Property<bool> m_signedradius {this, "SignedDriftRadius", true};
    Gaudi::Property<bool> m_calomat {this, "MuidMat", false};
    Gaudi::Property<bool> m_extmat {this, "ExtrapolatorMaterial", true};
    Gaudi::Property<bool> m_fillderivmatrix {this, "FillDerivativeMatrix", false};
    Gaudi::Property<bool> m_straightlineprop {this, "StraightLine", true};
    Gaudi::Property<bool> m_extensioncuts {this, "TRTExtensionCuts", true};
    Gaudi::Property<bool> m_sirecal {this, "RecalibrateSilicon", false};
    Gaudi::Property<bool> m_trtrecal {this, "RecalibrateTRT", false};
    Gaudi::Property<bool> m_kinkfinding {this, "KinkFinding", false};
    Gaudi::Property<bool> m_decomposesegments {this, "DecomposeSegments", true};
    Gaudi::Property<bool> m_getmaterialfromtrack {this, "GetMaterialFromTrack", true};
    Gaudi::Property<bool> m_domeastrackpar {this, "MeasuredTrackParameters", true};
    Gaudi::Property<bool> m_storemat {this, "StoreMaterialOnTrack", true};
    Gaudi::Property<bool> m_redoderivs {this, "RecalculateDerivatives", false};
    Gaudi::Property<bool> m_reintoutl {this, "ReintegrateOutliers", false};
    Gaudi::Property<bool> m_acceleration {this, "Acceleration", false};
    Gaudi::Property<bool> m_numderiv {this, "NumericalDerivs", false};
    Gaudi::Property<bool> m_fiteloss {this, "FitEnergyLoss", false};
    Gaudi::Property<bool> m_asymeloss {this, "AsymmetricEnergyLoss", true};
    Gaudi::Property<bool> m_useCaloTG {this, "UseCaloTG", false};
    Gaudi::Property<bool> m_rejectLargeNScat {this, "RejectLargeNScat", false};
    Gaudi::Property<bool> m_createSummary {this, "CreateTrackSummary", true};
    Gaudi::Property<bool> m_holeSearch {this, "DoHoleSearch", false};
 
    Gaudi::Property<double> m_outlcut {this, "OutlierCut", 5.0};
    Gaudi::Property<double> m_p {this, "Momentum", 0.0};
    Gaudi::Property<double> m_chi2cut {this, "TrackChi2PerNDFCut", 1.e15};
    Gaudi::Property<double> m_scalefactor {this, "TRTTubeHitCut", 2.5};
    Gaudi::Property<double> m_minphfcut {this, "MinPHFCut", 0.};
 
    Gaudi::Property<int> m_maxoutliers {this, "MaxOutliers", 10};
    Gaudi::Property<int> m_maxit {this, "MaxIterations", 30};
    Gaudi::Property<int> m_miniter {this, "MinimumIterations", 1};
    Gaudi::Property<int> m_fixbrem {this, "FixBrem", -1};
    Gaudi::Property<int> m_maxitPixelROT {this, "IterationsToRebuildPixelRots", 0};
 
    SG::ReadHandleKey<Trk::ClusterSplitProbabilityContainer>   m_clusterSplitProbContainer{this, "ClusterSplitProbabilityName", "",""};
 
 
    /*
     * This little volume defines the inner detector. Its exact size is set at
     * the time that the fitter object is created. We just make one and keep
     * it around to save us a few allocations.
     */
    Trk::Volume m_idVolume;
 
    /*
     * The following members are mutable. They keep track of the number of
     * fits that have returned with a certain status. Since this must be
     * shared across threads, we protect the array with a mutex, and we mark
     * these members as thread_safe for the ATLAS G++ plugin.
     */
    //mutable std::mutex m_fit_status_lock ATLAS_THREAD_SAFE;
    mutable std::array<std::atomic<unsigned int>, __S_MAX_VALUE> m_fit_status ATLAS_THREAD_SAFE = {};
  };
}
#endif