text: #ifndef MuonAnalysis_MuonAssociators_MatcherUsingTracksAlgorithm_h
  #define MuonAnalysis_MuonAssociators_MatcherUsingTracksAlgorithm_h
  //
  //
  
  /**
    \class    pat::MatcherUsingTracksAlgorithm MatcherUsingTracksAlgorithm.h "MuonAnalysis/MuonAssociators/interface/MatcherUsingTracksAlgorithm.h"
    \brief    Matcher of reconstructed objects to other reconstructed objects using the tracks inside them 
              
    \author   Giovanni Petrucciani
    \version  $Id: MatcherUsingTracksAlgorithm.h,v 1.7 2010/07/12 20:56:10 gpetrucc Exp $
  */
  
  #include "DataFormats/Candidate/interface/CandidateFwd.h"
  #include "DataFormats/TrackReco/interface/TrackFwd.h"
  
  #include "FWCore/Framework/interface/ESHandle.h"
  #include "FWCore/Framework/interface/EventSetup.h"
  #include "FWCore/Framework/interface/ConsumesCollector.h"
  #include "FWCore/Utilities/interface/InputTag.h"
  #include "FWCore/ParameterSet/interface/ParameterSet.h"
  
  #include "Geometry/CommonDetUnit/interface/GlobalTrackingGeometry.h"
  #include "MagneticField/Engine/interface/MagneticField.h"
  #include "TrackingTools/GeomPropagators/interface/Propagator.h"
  #include "TrackingTools/TrajectoryState/interface/TrajectoryStateOnSurface.h"
  #include "TrackingTools/TrajectoryState/interface/TrajectoryStateClosestToPoint.h"
  #include "CommonTools/Utils/interface/StringCutObjectSelector.h"
  
  class IdealMagneticFieldRecord;
  class TrackingComponentsRecord;
  class GlobalTrackingGeometryRecord;
  
  class MatcherUsingTracksAlgorithm {
  public:
    explicit MatcherUsingTracksAlgorithm(const edm::ParameterSet &iConfig, edm::ConsumesCollector);
    virtual ~MatcherUsingTracksAlgorithm() {}
  
    /// Call this method at the beginning of each run, to initialize geometry, magnetic field and propagators
    void init(const edm::EventSetup &iSetup);
  
    /// Try to match one track to another one. Return true if succeeded.
    /// For matches not by ref, it will update deltaR, deltaLocalPos and deltaPtRel if the match suceeded
    bool match(const reco::Candidate &c1,
               const reco::Candidate &c2,
               float &deltaR,
               float &deltaEta,
               float &deltaPhi,
               float &deltaLocalPos,
               float &deltaPtRel,
               float &chi2) const;
  
    /// Find the best match to another candidate, and return its index in the vector
    /// For matches not by ref, it will update deltaR, deltaLocalPos and deltaPtRel if the match suceeded
    /// Returns -1 if the match fails
    int match(const reco::Candidate &tk,
              const edm::View<reco::Candidate> &c2s,
              float &deltaR,
              float &deltaEta,
              float &deltaPhi,
              float &deltaLocalPos,
              float &deltaPtRel,
              float &chi2) const;
  
    /// Return 'true' if the matcher will produce meaningful deltaR, deltaLocalPos, deltaPtRel values
    bool hasMetrics() const { return algo_ != ByTrackRef; }
  
    /// Return 'true' if the matcher will produce also chi2
    bool hasChi2() const { return useChi2_; }
  
    /// Compute the chi2 of two free trajectory states, in the curvilinear frame (q/p, theta, phi, dxy, dsz)
    /// At least one must have errors
    /// diagonalOnly: don't use off-diagonal terms of covariance matrix
    /// useVertex   : use dxy, dsz in the chi2 (if false, use only q/p, theta, phi)
    /// useFirstMomentum : use the 'start' state momentum to compute dxy, dsx (if false, use 'other')
    static double getChi2(const FreeTrajectoryState &start,
                          const FreeTrajectoryState &other,
                          bool diagonalOnly,
                          bool useVertex,
                          bool useFirstMomentum);
  
    /// Compute the chi2 of one free trajectory state and a TrajectoryStateClosestToPoint closest to it, in the perigee frame
    /// At least one must have errors
    /// diagonalOnly: don't use off-diagonal terms of covariance matrix
    /// useVertex   : use dxy, dsz in the chi2 (if false, use only q/p, theta, phi)
    static double getChi2(const FreeTrajectoryState &start,
                          const TrajectoryStateClosestToPoint &other,
                          bool diagonalOnly,
                          bool useVertex);
  
    /// Compute the chi2 of two free trajectory states, in the local frame (q/p, dx, dy, dxdz, dydz)
    /// At least one must have errors
    /// diagonalOnly: don't use off-diagonal terms of covariance matrix
    /// useVertex   : use dx, dy in the chi2 (if false, use only direction and q/p)
    static double getChi2(const TrajectoryStateOnSurface &start,
                          const TrajectoryStateOnSurface &other,
                          bool diagonalOnly,
                          bool usePosition);
  
    /// Possibly crop the 3x3 part of the matrix or remove off-diagonal terms, then invert.
    static void cropAndInvert(AlgebraicSymMatrix55 &cov, bool diagonalOnly, bool top3by3only);
  
  private:
    enum AlgoType {
      ByTrackRef,
      ByDirectComparison,
      ByPropagatingSrc,
      ByPropagatingMatched,
      ByPropagatingSrcTSCP,
      ByPropagatingMatchedTSCP
    };  // propagate closest to point
    enum WhichTrack { None, TrackerTk, MuonTk, GlobalTk };
    enum WhichState { AtVertex, Innermost, Outermost };
  
    AlgoType algo_;
    WhichTrack whichTrack1_, whichTrack2_;
    WhichState whichState1_, whichState2_;
  
    // Preselection cuts on both sides
    StringCutObjectSelector<reco::Candidate, true> srcCut_, matchedCut_;
  
    // Matching cuts
    float maxLocalPosDiff_;
    float maxGlobalMomDeltaR_;
    float maxGlobalMomDeltaEta_;
    float maxGlobalMomDeltaPhi_;
    float maxGlobalDPtRel_;
    float maxChi2_;
    bool requireSameCharge_;
    bool useChi2_, chi2UseVertex_, chi2DiagonalOnly_, chi2FirstMomentum_;
    enum SortBy { LocalPosDiff, GlobalMomDeltaR, GlobalMomDeltaEta, GlobalMomDeltaPhi, GlobalDPtRel, Chi2 };
    SortBy sortBy_;
  
    //- Tools
    edm::ESHandle<MagneticField> magfield_;
    edm::ESHandle<Propagator> propagator_;
    edm::ESHandle<GlobalTrackingGeometry> geometry_;
  
    edm::ESGetToken<MagneticField, IdealMagneticFieldRecord> magfieldToken_;
    edm::ESGetToken<Propagator, TrackingComponentsRecord> propagatorToken_;
    edm::ESGetToken<GlobalTrackingGeometry, GlobalTrackingGeometryRecord> geometryToken_;
  
    /// Get track reference out of a Candidate (via dynamic_cast to reco::RecoCandidate)
    reco::TrackRef getTrack(const reco::Candidate &reco, WhichTrack which) const;
  
    /// Starting state for the propagation
    FreeTrajectoryState startingState(const reco::Candidate &reco, WhichTrack whichTrack, WhichState whichState) const;
  
    /// End state for the propagation
    TrajectoryStateOnSurface targetState(const reco::Candidate &reco, WhichTrack whichTrack, WhichState whichState) const;
  
    /// Propagate and match. return true if current pair is the new best match (in that case, update also deltaR and deltaLocalPos)
    /// Uses TrajectoryStateClosestToPointBuilder
    bool matchWithPropagation(const FreeTrajectoryState &start,
                              const FreeTrajectoryState &target,
                              float &lastDeltaR,
                              float &lastDeltaEta,
                              float &lastDeltaPhi,
                              float &lastDeltaLocalPos,
                              float &lastGlobalDPtRel,
                              float &lastChi2) const;
  
    /// Propagate and match. return true if current pair is the new best match (in that case, update also deltaR and deltaLocalPos)
    /// Uses standard propagator to reach target's surface
    bool matchWithPropagation(const FreeTrajectoryState &start,
                              const TrajectoryStateOnSurface &target,
                              float &lastDeltaR,
                              float &lastDeltaEta,
                              float &lastDeltaPhi,
                              float &lastDeltaLocalPos,
                              float &lastGlobalDPtRel,
                              float &lastChi2) const;
  
    /// Compare directly two states. return true if current pair is the new best match (in that case, update also deltaR and deltaLocalPos)
    bool matchByDirectComparison(const FreeTrajectoryState &start,
                                 const FreeTrajectoryState &other,
                                 float &lastDeltaR,
                                 float &lastDeltaEta,
                                 float &lastDeltaPhi,
                                 float &lastDeltaLocalPos,
                                 float &lastGlobalDPtRel,
                                 float &lastChi2) const;
  
    /// Parse some configuration
    void getConf(const edm::ParameterSet &iConfig,
                 const std::string &whatFor,
                 WhichTrack &whichTrack,
                 WhichState &whichState);
  };
  
  #endif
  