VrtSecInclusive.h

Go to the documentation of this file.

/*
  Copyright (C) 2002-2024 CERN for the benefit of the ATLAS collaboration
*/
 
// VKalVrt.h
//
#ifndef _VrtSecInclusive_VrtSecInclusive_H
#define _VrtSecInclusive_VrtSecInclusive_H
 
 
#include "VrtSecInclusive/Constants.h"
 
#include "AthenaBaseComps/AthAlgorithm.h"
 
// Gaudi includes
#include "GaudiKernel/ToolHandle.h"
#include "GaudiKernel/ITHistSvc.h"
//
#include "TrkVKalVrtFitter/TrkVKalVrtFitter.h"
 
// for truth
#include "GeneratorObjects/McEventCollection.h"
 
#include "TrkToolInterfaces/ITruthToTrack.h"
#include "ITrackToVertex/ITrackToVertex.h"
#include "TrkVertexFitterInterfaces/ITrackToVertexIPEstimator.h"
#include "TrkExInterfaces/IExtrapolator.h"
#include "TrkExInterfaces/IPropagator.h"
#include "TrkSurfaces/CylinderSurface.h"
#include "TrkDetDescrInterfaces/IVertexMapper.h"
#include "GaudiKernel/ServiceHandle.h"
#include "InDetConditionsSummaryService/IInDetConditionsTool.h"
#include "InDetIdentifier/PixelID.h"
#include "InDetIdentifier/SCT_ID.h"
 
// xAOD Classes
#include "xAODEventInfo/EventInfo.h"
#include "xAODTracking/TrackParticleContainer.h"
#include "xAODTracking/TrackParticleAuxContainer.h"
#include "xAODTracking/TrackParticle.h"
#include "xAODTracking/VertexContainer.h"
#include "xAODTracking/VertexAuxContainer.h"
#include "xAODTruth/TruthParticleContainer.h"
#include "xAODTruth/TruthEventContainer.h"
#include "xAODTruth/TruthVertexContainer.h"
#include "xAODMuon/MuonContainer.h"
#include "xAODMuon/Muon.h"
#include "xAODEgamma/ElectronContainer.h"
#include "xAODEgamma/Electron.h"
 
// Normal STL and physical vectors
#include <vector>
#include <deque>
#include <functional>
 
 
class TH1;
 
namespace Trk {
  class ITruthToTrack;
  class ITrackToVertex;
  class ITrackToVertexIPEstimator;
  class IExtrapolator;
  class IVertexMapper;
  struct MappedVertex;
}
 
namespace VKalVrtAthena {
 
  namespace GeoModel { enum GeoModel { Run1=1, Run2=2 }; }
 
  class IntersectionPos;
  class IntersectionPos_barrel;
  class IntersectionPos_endcap;
 
  class NtupleVars;
}
 
 
namespace VKalVrtAthena {
 
  class VrtSecInclusive : public AthAlgorithm {
  public:
    VrtSecInclusive(const std::string& name, ISvcLocator* pSvcLocator);
 
    ~VrtSecInclusive();
 
    virtual StatusCode initialize();
    virtual StatusCode finalize();
    virtual StatusCode execute();
    virtual StatusCode initEvent();
 
  private:
 
    //
    //  Member Variables
    //
 
    struct JobProperties {
      // JO: GeoModel
      int    geoModel;
 
      std::string          TrackLocation;
      std::string          MuonLocation;
      std::string          ElectronLocation;
      std::string          PrimVrtLocation;
      std::string          truthParticleContainerName;
      std::string          mcEventContainerName;
      std::string          augVerString;
      std::string          truthParticleFilter;
 
      std::string all2trksVerticesContainerName;
      std::string secondaryVerticesContainerName;
 
      // common feature flags
      bool   doTruth;
      bool   FillHist;
      bool   FillNtuple;
      bool   FillIntermediateVertices;
      bool   doIntersectionPos;
      bool   doMapToLocal;
      bool   extrapPV; //extrapolate reco and prim truth trks to PV (for testing only)
 
      bool   passThroughTrackSelection;
 
      bool   doFastMode; // flag for running in rapid finder mode instead of using graph
 
      // track selection conditions
      unsigned int SelTrkMaxCutoff;
      bool  SAloneTRT;
 
      /* impact parameters */
      bool   do_PVvetoCut;
      bool   do_d0Cut;
      bool   do_z0Cut;
      bool   do_d0errCut;
      bool   do_z0errCut;
      bool   do_d0signifCut;
      bool   do_z0signifCut;
 
      bool   ImpactWrtBL;
      double d0TrkPVDstMinCut;
      double d0TrkPVDstMaxCut;
      double d0TrkPVSignifCut;
      double z0TrkPVDstMinCut;
      double z0TrkPVDstMaxCut;
      double z0TrkPVSignifCut;
      double d0TrkErrorCut;
      double z0TrkErrorCut;
      double twoTrkVtxFormingD0Cut;
 
      /* pT anc chi2 */
      double TrkChi2Cut;
      double TrkPtCut;
 
      /* hit requirements */
      bool doTRTPixCut; // Kazuki
      int  CutSctHits;
      int  CutPixelHits;
      int  CutSiHits;
      int  CutBLayHits;
      int  CutSharedHits;
      int  CutTRTHits; // Kazuki
      int  CutTightSCTHits;
      int  CutTightTRTHits;
 
      /* track extrpolator; 1==VKalGetImpact, 2==m_trackToVertexTool*/
      int trkExtrapolator;
 
      // Vertex reconstruction
      bool   doPVcompatibilityCut;
      bool   doTightPVcompatibilityCut;
      bool   removeFakeVrt;
      bool   removeFakeVrtLate;
      bool   doReassembleVertices;
      bool   doMergeByShuffling;
      bool   doSuggestedRefitOnMerging;
      bool   doMagnetMerging;
      bool   doWildMerging;
      bool   doMergeFinalVerticesDistance; // Kazuki
      bool   doAssociateNonSelectedTracks;
      bool   doFinalImproveChi2;
      double pvCompatibilityCut;
      double SelVrtChi2Cut;
      double VertexMergeFinalDistCut; // Kazuki
      double VertexMergeFinalDistScaling;
      double VertexMergeCut;
      double TrackDetachCut;
 
      bool   doTwoTrSoftBtag;
      double twoTrVrtAngleCut;
      double twoTrVrtMinDistFromPV;
 
      // When truncateWrkVertices set to true, maxWrkVertices is the maximum
      // number of potential vertices to process, used to truncate the list
      // in rare caes where several thousands are found, to avoid algorithm
      // timeout issues
      bool   truncateWrkVertices;
      size_t maxWrkVertices;
 
      double associateMinDistanceToPV;
      double associateMaxD0Signif;
      double associateMaxZ0Signif;
      double associatePtCut;
      double associateChi2Cut;
 
      double reassembleMaxImpactParameterD0;
      double reassembleMaxImpactParameterZ0;
      double mergeByShufflingMaxSignificance;
      double mergeByShufflingAllowance;
 
      double improveChi2ProbThreshold;
 
      // vertexing using muons (test implementation)
      bool doSelectTracksFromMuons;
      bool doRemoveCaloTaggedMuons;
      bool doSelectTracksFromElectrons;
      bool doSelectIDAndGSFTracks;
      bool doRemoveNonLeptonVertices;
 
      // When doSelectTracksWithLRTCuts is set to true, the addtional track cuts
      // be applied to the selected tracks to reduce the number of fake tracks in
      // the selected track collected. These cuts are inspired by the improvments that
      // were implmented for LRT Run 3.
      bool doSelectTracksWithLRTCuts ;
 
      // Additional dressing option
      bool doAugmentDVimpactParametersToMuons;     // potentially useful for DV + muon search
      bool doAugmentDVimpactParametersToElectrons; // potentially useful for analyses involving electrons
 
      // MC truth
      double               mcTrkResolution;
      double               TruthTrkLen;
 
    };
 
    struct JobProperties m_jp;
 
    // Indicates give-up modes during vertexing
    // 0 if no errors occured
    // 1 if too few selected tracks
    // 2 if too many selected tracks
    // 3 if wrkVertices container is truncated
    // 4 if any uncaught exception is raised at the top level of VSI execute()
    int m_vertexingStatus = 0;
 
    // xAOD Accessors
    const xAOD::VertexContainer*  m_primaryVertices;
    const xAOD::Vertex*           m_thePV;
    std::unique_ptr<std::vector<const xAOD::TrackParticle*> > m_selectedTracks;
    std::unique_ptr<std::vector<const xAOD::TrackParticle*> > m_associatedTracks;
    std::unique_ptr<std::vector<const xAOD::TrackParticle*> > m_leptonicTracks;
    std::vector<double>  m_BeamPosition;
 
    //
    //  Athena JobOption Properties
    //
 
    ToolHandle <Trk::ITrkVKalVrtFitter>        m_fitSvc;       // VKalVrtFitter tool
    ToolHandle <Trk::ITruthToTrack>            m_truthToTrack; // tool to create trkParam from genPart
 
    ToolHandle< Reco::ITrackToVertex >         m_trackToVertexTool;
    ToolHandle<Trk::ITrackToVertexIPEstimator> m_trackToVertexIPEstimatorTool;
    ToolHandle<Trk::IExtrapolator>             m_extrapolator;
    ToolHandle<Trk::IVertexMapper>             m_vertexMapper;
 
    ToolHandle<IInDetConditionsTool> m_pixelCondSummaryTool{this, "PixelConditionsSummaryTool", "PixelConditionsSummaryTool", "Tool to retrieve Pixel Conditions summary"};
    ToolHandle<IInDetConditionsTool> m_sctCondSummaryTool{this, "InDetSCT_ConditionsSummaryTool", "SCT_ConditionsSummaryTool/InDetSCT_ConditionsSummaryTool", "Tool to retrieve SCT conditions summary"};
 
    const AtlasDetectorID* m_atlasId = nullptr;
    const PixelID* m_pixelId = nullptr;
    const SCT_ID*  m_sctId = nullptr;
 
    std::string m_checkPatternStrategy;
    using PatternStrategyFunc = bool (VrtSecInclusive::*) ( const xAOD::TrackParticle *trk, const Amg::Vector3D& vertex );
    std::map<std::string, PatternStrategyFunc> m_patternStrategyFuncs;
 
    // AuxElement decorators
    std::unique_ptr< SG::AuxElement::Decorator< char > > m_decor_isSelected;
    std::unique_ptr< SG::AuxElement::Decorator< char > > m_decor_isAssociated;
    std::unique_ptr< SG::AuxElement::Decorator< char > > m_decor_is_svtrk_final;
    std::map< unsigned, SG::AuxElement::Decorator<float> > m_trkDecors;
 
    SG::ReadHandleKey<xAOD::EventInfo> m_eventInfoKey{this,"EventInfoKey", "EventInfo", "EventInfo name"};
    SG::WriteDecorHandleKey<xAOD::EventInfo> m_vertexingStatusKey;
 
    using IPDecoratorType = SG::AuxElement::Decorator< std::vector< std::vector<float> > >;
    std::unique_ptr< IPDecoratorType > m_decor_d0_wrtSVs;
    std::unique_ptr< IPDecoratorType > m_decor_z0_wrtSVs;
    std::unique_ptr< IPDecoratorType > m_decor_pt_wrtSVs;
    std::unique_ptr< IPDecoratorType > m_decor_eta_wrtSVs;
    std::unique_ptr< IPDecoratorType > m_decor_phi_wrtSVs;
    std::unique_ptr< IPDecoratorType > m_decor_d0err_wrtSVs;
    std::unique_ptr< IPDecoratorType > m_decor_z0err_wrtSVs;
 
    using VertexELType = SG::AuxElement::Decorator< std::vector<ElementLink< xAOD::VertexContainer > > >;
    std::unique_ptr< VertexELType > m_decor_svLink;
 
    //
    // define ntuple variables here
    //
 
    // The standard AANT, CollectionTree, is bare bones
    TTree      *m_tree_Vert;
    std::unique_ptr<NtupleVars> m_ntupleVars;
 
    // Histograms for stats
    std::map<std::string, TH1*> m_hists;
 
 
    //
    // Private member functions
    //
 
    // for event info to new ntuple (used to go by default in CollectionTree)
    void declareProperties();
 
    StatusCode addEventInfo();
    StatusCode setupNtupleVariables();
    StatusCode setupNtuple();
    StatusCode clearNtupleVariables();
    StatusCode deleteNtupleVariables();
    StatusCode processPrimaryVertices();
    StatusCode fillAANT_SelectedBaseTracks();
    StatusCode fillAANT_SecondaryVertices( xAOD::VertexContainer* );
 
    //
    struct WrkVrt { 
      bool isGood = false;                            
      std::deque<long int> selectedTrackIndices;      
      std::deque<long int> associatedTrackIndices;    
      Amg::Vector3D        vertex;                    
      TLorentzVector       vertexMom;                 
      std::vector<double>  vertexCov;                 
      double               Chi2 = 0;                  
      double               Chi2_core = 0;             
      std::vector<double>  Chi2PerTrk;                
      long int             Charge = 0;                
      std::vector< std::vector<double> > TrkAtVrt;    
      unsigned long        closestWrkVrtIndex = 0;    
      double               closestWrkVrtValue = 0;    
 
      inline double ndof() const { return 2.0*( selectedTrackIndices.size() + associatedTrackIndices.size() ) - 3.0; }
      inline double ndof_core() const { return 2.0*( selectedTrackIndices.size() ) - 3.0; }
      inline unsigned nTracksTotal() const { return selectedTrackIndices.size() + associatedTrackIndices.size(); }
      inline double fitQuality() const { return Chi2 / ndof(); }
    };
 
 
    using Detector = int;
    using Bec      = int;
    using Layer    = int;
    using Flag     = int;
    using ExtrapolatedPoint   = std::tuple<const TVector3, Detector, Bec, Layer, Flag>;
    using ExtrapolatedPattern = std::vector< ExtrapolatedPoint >;
    using PatternBank         = std::map<const xAOD::TrackParticle*, std::pair< std::unique_ptr<ExtrapolatedPattern>, std::unique_ptr<ExtrapolatedPattern> > >;
 
    PatternBank m_extrapolatedPatternBank;
 
    std::vector< std::pair<int, int> > m_incomp;
 
    // the map used by printWrkSet
    std::map<const xAOD::TruthVertex*, bool> m_matchMap;
 
 
    void selectTrack( const xAOD::TrackParticle* );
    StatusCode selectTracksInDet();
    StatusCode selectTracksFromMuons();
    StatusCode selectTracksFromElectrons();
    StatusCode selectInDetAndGSFTracks();
 
    using TrackSelectionAlg = StatusCode (VrtSecInclusive::*)();
    std::vector<TrackSelectionAlg> m_trackSelectionAlgs;
 
    using CutFunc = bool (VrtSecInclusive::*) ( const xAOD::TrackParticle* ) const;
    std::vector<CutFunc> m_trackSelectionFuncs;
 
    bool selectTrack_notPVassociated ( const xAOD::TrackParticle* ) const;
    bool selectTrack_pTCut           ( const xAOD::TrackParticle* ) const;
    bool selectTrack_chi2Cut         ( const xAOD::TrackParticle* ) const;
    bool selectTrack_hitPattern      ( const xAOD::TrackParticle* ) const;
    bool selectTrack_hitPatternTight ( const xAOD::TrackParticle* ) const;
    bool selectTrack_d0Cut           ( const xAOD::TrackParticle* ) const;
    bool selectTrack_z0Cut           ( const xAOD::TrackParticle* ) const;
    bool selectTrack_d0errCut        ( const xAOD::TrackParticle* ) const;
    bool selectTrack_z0errCut        ( const xAOD::TrackParticle* ) const;
    static bool selectTrack_d0signifCut     ( const xAOD::TrackParticle* ) ;
    static bool selectTrack_z0signifCut     ( const xAOD::TrackParticle* ) ;
    bool selectTrack_LRTR3Cut        ( const xAOD::TrackParticle* ) const;
 
    StatusCode extractIncompatibleTrackPairs( std::vector<WrkVrt>* );
    StatusCode findNtrackVertices(std::vector<WrkVrt>* );
    StatusCode rearrangeTracks( std::vector<WrkVrt>* );
 
    StatusCode reassembleVertices( std::vector<WrkVrt>* );
 
    StatusCode mergeByShuffling( std::vector<WrkVrt>* );
 
    StatusCode mergeFinalVertices( std::vector<WrkVrt>* ); // Kazuki
 
    StatusCode associateNonSelectedTracks( std::vector<WrkVrt>* );
 
    StatusCode refitAndSelectGoodQualityVertices( std::vector<WrkVrt>* );
 
    bool getSVImpactParameters(const xAOD::TrackParticle* trk, const Amg::Vector3D& vertex, std::vector<double>& impactParameters, std::vector<double>& impactParErrors);
 
    enum TrkParameter    { k_d0=0, k_z0=1, k_theta=2, k_phi=3, k_qOverP=4 ,k_nTP=5 };
    enum TrkParameterUnc { k_d0d0=0, k_z0z0=1, k_nTPU=2 };
 
    using vertexingAlg = StatusCode (VrtSecInclusive::*)( std::vector<WrkVrt>* );
    std::vector< std::pair<std::string, vertexingAlg> > m_vertexingAlgorithms;
    unsigned m_vertexingAlgorithmStep = 0U;
 
 
    //
    // Supporting utility functions
 
    void printWrkSet(const std::vector<WrkVrt> *WrkVrtSet, const std::string& name);
 
    StatusCode refitVertex( WrkVrt& );
    StatusCode refitVertex( WrkVrt&, Trk::IVKalState& istate );
 
    StatusCode refitVertexWithSuggestion( WrkVrt&, const Amg::Vector3D& );
    StatusCode refitVertexWithSuggestion( WrkVrt&, const Amg::Vector3D&, Trk::IVKalState& istate );
 
    double improveVertexChi2( WrkVrt& );
 
    static void removeTrackFromVertex(std::vector<WrkVrt>*,
                               std::vector< std::deque<long int> > *,
             const long int & ,const long int & );
 
    StatusCode disassembleVertex(std::vector<WrkVrt> *, const unsigned& vertexIndex );
 
    void trackClassification(std::vector< WrkVrt >* , std::map< long int, std::vector<long int> >& );
 
    double findWorstChi2ofMaximallySharedTrack(std::vector<WrkVrt>*, std::map< long int, std::vector<long int> >&, long int & ,long int & );
 
    static size_t nTrkCommon( std::vector<WrkVrt> *WrkVrtSet, const std::pair<unsigned, unsigned>& pairIndex ) ;
 
    double significanceBetweenVertices( const WrkVrt&, const WrkVrt& ) const;
 
    double distanceBetweenVertices( const WrkVrt&, const WrkVrt& ) const;
 
    using AlgForVerticesPair = double (VrtSecInclusive::*)( const WrkVrt&, const WrkVrt& ) const;
 
    double findMinVerticesPair( std::vector<WrkVrt>*, std::pair<unsigned, unsigned>&, const AlgForVerticesPair& );
    
    static double findMinVerticesNextPair( std::vector<WrkVrt>*, std::pair<unsigned, unsigned>& );
 
    StatusCode mergeVertices( WrkVrt& destination, WrkVrt& source );
 
    enum mergeStep { RECONSTRUCT_NTRK, REASSEMBLE, SHUFFLE1, SHUFFLE2, SHUFFLE3, FINAL };
 
    typedef struct track_summary_properties {
      uint8_t numIBLHits;
      uint8_t numBLayerHits;
      uint8_t numPixelLayer1_Hits;
      uint8_t numPixelLayer2_Hits;
      uint8_t numPixelDisk0_Hits;
      uint8_t numPixelDisk1_Hits;
      uint8_t numPixelDisk2_Hits;
      uint8_t numPixelHits;
      uint8_t numSctBarrelLayer0_Hits;
      uint8_t numSctBarrelLayer1_Hits;
      uint8_t numSctBarrelLayer2_Hits;
      uint8_t numSctBarrelLayer3_Hits;
      uint8_t numSctEC0_Hits;
      uint8_t numSctEC1_Hits;
      uint8_t numSctEC2_Hits;
      uint8_t numSctEC3_Hits;
      uint8_t numSctEC4_Hits;
      uint8_t numSctHits;
      uint8_t numTrtHits;
    } track_summary;
 
    static void fillTrackSummary( track_summary& summary, const xAOD::TrackParticle *trk );
 
    ExtrapolatedPattern* extrapolatedPattern( const xAOD::TrackParticle*, enum Trk::PropDirection );
 
    bool patternCheck    ( const uint32_t& pattern, const Amg::Vector3D& vertex );
    static bool patternCheckRun1( const uint32_t& pattern, const Amg::Vector3D& vertex );
    static bool patternCheckRun2( const uint32_t& pattern, const Amg::Vector3D& vertex );
 
    bool patternCheckOuterOnly    ( const uint32_t& pattern, const Amg::Vector3D& vertex );
    static bool patternCheckRun1OuterOnly( const uint32_t& pattern, const Amg::Vector3D& vertex );
    static bool patternCheckRun2OuterOnly( const uint32_t& pattern, const Amg::Vector3D& vertex );
 
    bool checkTrackHitPatternToVertex( const xAOD::TrackParticle *trk, const Amg::Vector3D& vertex );
 
    bool checkTrackHitPatternToVertexOuterOnly( const xAOD::TrackParticle *trk, const Amg::Vector3D& vertex );
 
    bool checkTrackHitPatternToVertexByExtrapolation( const xAOD::TrackParticle *trk, const Amg::Vector3D& vertex );
 
    bool checkTrackHitPatternToVertexByExtrapolationAssist( const xAOD::TrackParticle *trk, const Amg::Vector3D& vertex );
 
    bool passedFakeReject( const Amg::Vector3D& FitVertex, const xAOD::TrackParticle *itrk, const xAOD::TrackParticle *jtrk );
 
    void removeInconsistentTracks( WrkVrt& );
 
    template<class Track> void getIntersection(Track *trk, std::vector<IntersectionPos*>& layers, const Trk::Perigee* per);
    template<class Track> void setIntersection(Track *trk, IntersectionPos *bec, const Trk::Perigee* per);
 
    StatusCode monitorVertexingAlgorithmStep( std::vector<WrkVrt>*, const std::string& name, bool final = false );
 
    //
    // Truth Information Algorithms Member Functions
    //
    //
 
    static const xAOD::TruthParticle *getTrkGenParticle(const xAOD::TrackParticle*) ;
 
    StatusCode categorizeVertexTruthTopology( xAOD::Vertex *vertex );
 
    void dumpTruthInformation();
 
    std::vector<const xAOD::TruthVertex*> m_tracingTruthVertices;
 
    //
    // Additional augmentation
    //
    //
 
    template<class LeptonFlavor>
    StatusCode augmentDVimpactParametersToLeptons( const std::string& containerName );
 
  };
 
} // end of namespace bracket
 
 
// This header file contains the definition of member templates
#include "details/Utilities.h"
 
 
#endif /* _VrtSecInclusive_VrtSecInclusive_H */