DuplicateTrackMerger.cc

Go to the documentation of this file.

#include "FWCore/Framework/interface/stream/EDProducer.h"

#include "FWCore/ParameterSet/interface/ParameterSet.h"
#include "FWCore/ParameterSet/interface/ParameterSetDescription.h"
#include "FWCore/ParameterSet/interface/ConfigurationDescriptions.h"
#include "FWCore/MessageLogger/interface/MessageLogger.h"


#include "DataFormats/TrackReco/interface/TrackFwd.h"
#include "DataFormats/TrackReco/interface/Track.h"
#include "FWCore/MessageLogger/interface/MessageLogger.h"

#include "DataFormats/TrackCandidate/interface/TrackCandidate.h"

#include "TrackingTools/PatternTools/interface/TSCPBuilderNoMaterial.h"
#include "TrackingTools/GeomPropagators/interface/Propagator.h"
#include "TrackingTools/KalmanUpdators/interface/Chi2MeasurementEstimatorBase.h"
#include "TrackMerger.h"

#include <vector>
#include <algorithm>
#include <string>
#include <iostream>
#include <atomic>

#include "CondFormats/EgammaObjects/interface/GBRForest.h"

using namespace reco;
namespace {
  class  DuplicateTrackMerger final : public edm::stream::EDProducer<> {
  public:
    explicit DuplicateTrackMerger(const edm::ParameterSet& iPara);
    virtual ~DuplicateTrackMerger();
    
    using CandidateToDuplicate = std::vector<std::pair<int, int>>;
     
    static void fillDescriptions(edm::ConfigurationDescriptions& descriptions);

     
    private:
      void produce( edm::Event &, const edm::EventSetup &) override;
      
      bool checkForDisjointTracks(const reco::Track *t1, const reco::Track *t2, TSCPBuilderNoMaterial& tscpBuilder) const;
      bool checkForOverlappingTracks(const reco::Track *t1, const reco::Track *t2, const float cosT) const;

    private:
      const GBRForest* forest_;
      
           
      std::string dbFileName_;
      bool useForestFromDB_;
      std::string forestLabel_;
      
      std::string propagatorName_;
      std::string chi2EstimatorName_;

      edm::EDGetTokenT<reco::TrackCollection> trackSource_;
      double minDeltaR3d2_;
      double minBDTG_;
      double minpT2_;
      double minP_;
      float maxDCA2_;
      float maxDPhi_;
      float maxDLambda_;
      float maxDdxy_;
      float maxDdsz_;
      float maxDQoP_;
      unsigned int overlapCheckMaxHits_;
      unsigned int overlapCheckMaxMissingLayers_;
      float overlapCheckMinCosT_;

      const MagneticField *magfield_;
      const TrackerTopology *ttopo_;
      const TrackerGeometry *geom_;
      const Propagator *propagator_;
      const Chi2MeasurementEstimatorBase *chi2Estimator_;

      TrackMerger merger_;

#ifdef EDM_ML_DEBUG
      bool debug_;
#endif
    };
  }
    
#include "FWCore/Framework/interface/Event.h"
#include "TrackingTools/TrajectoryState/interface/FreeTrajectoryState.h"
#include "TrackingTools/TrajectoryState/interface/TrajectoryStateOnSurface.h"
#include "TrackingTools/TrajectoryState/interface/TrajectoryStateTransform.h"
#include "CondFormats/DataRecord/interface/GBRWrapperRcd.h"
#include "TrackingTools/Records/interface/TrackingComponentsRecord.h"
#include "FWCore/Framework/interface/EventSetup.h"
#include "FWCore/Framework/interface/ESHandle.h" 
#include "TFile.h"

#include "DataFormats/SiPixelDetId/interface/PixelSubdetector.h"
#include "DataFormats/TrackerRecHit2D/interface/BaseTrackerRecHit.h"

#include "DuplicateTrackType.h"

namespace {


void
DuplicateTrackMerger::fillDescriptions(edm::ConfigurationDescriptions& descriptions)
{
     edm::ParameterSetDescription desc;
     desc.add<edm::InputTag>("source",edm::InputTag());
     desc.add<double>("minDeltaR3d",-4.0);
     desc.add<double>("minBDTG",-0.1);
     desc.add<double>("minpT",0.2);
     desc.add<double>("minP",0.4);
     desc.add<double>("maxDCA",30.0);
     desc.add<double>("maxDPhi",0.30);
     desc.add<double>("maxDLambda",0.30);
     desc.add<double>("maxDdsz",10.0);
     desc.add<double>("maxDdxy",10.0);
     desc.add<double>("maxDQoP",0.25);
     desc.add<unsigned>("overlapCheckMaxHits", 4);
     desc.add<unsigned>("overlapCheckMaxMissingLayers", 1);
     desc.add<double>("overlapCheckMinCosT", 0.99);
     desc.add<std::string>("forestLabel","MVADuplicate");
     desc.add<std::string>("GBRForestFileName","");
     desc.add<bool>("useInnermostState",true);
     desc.add<std::string>("ttrhBuilderName","WithAngleAndTemplate");
     desc.add<std::string>("propagatorName", "PropagatorWithMaterial");
     desc.add<std::string>("chi2EstimatorName", "DuplicateTrackMergerChi2Est");
     descriptions.add("DuplicateTrackMerger", desc);
}

  
DuplicateTrackMerger::DuplicateTrackMerger(const edm::ParameterSet& iPara) : forest_(nullptr), merger_(iPara)
{

  trackSource_ = consumes<reco::TrackCollection>(iPara.getParameter<edm::InputTag>("source"));
  minDeltaR3d2_ = iPara.getParameter<double>("minDeltaR3d"); minDeltaR3d2_*=std::abs(minDeltaR3d2_);
  minBDTG_ = iPara.getParameter<double>("minBDTG");
  minpT2_ = iPara.getParameter<double>("minpT"); minpT2_ *= minpT2_;
  minP_ = iPara.getParameter<double>("minP");
  maxDCA2_ = iPara.getParameter<double>("maxDCA"); maxDCA2_*=maxDCA2_;
  maxDPhi_ = iPara.getParameter<double>("maxDPhi");
  maxDLambda_ = iPara.getParameter<double>("maxDLambda");
  maxDdsz_ = iPara.getParameter<double>("maxDdsz");
  maxDdxy_ = iPara.getParameter<double>("maxDdxy");
  maxDQoP_ = iPara.getParameter<double>("maxDQoP");
  overlapCheckMaxHits_ = iPara.getParameter<unsigned>("overlapCheckMaxHits");
  overlapCheckMaxMissingLayers_ = iPara.getParameter<unsigned>("overlapCheckMaxMissingLayers");
  overlapCheckMinCosT_ = iPara.getParameter<double>("overlapCheckMinCosT");

  produces<std::vector<TrackCandidate> >("candidates");
  produces<CandidateToDuplicate>("candidateMap");

  forestLabel_ = iPara.getParameter<std::string>("forestLabel");

  dbFileName_ = iPara.getParameter<std::string>("GBRForestFileName");
  useForestFromDB_ = dbFileName_.empty();

  propagatorName_ = iPara.getParameter<std::string>("propagatorName");
  chi2EstimatorName_ = iPara.getParameter<std::string>("chi2EstimatorName");

  /*
  tmvaReader_ = new TMVA::Reader("!Color:Silent");
  tmvaReader_->AddVariable("ddsz",&tmva_ddsz_);
  tmvaReader_->AddVariable("ddxy",&tmva_ddxy_);
  tmvaReader_->AddVariable("dphi",&tmva_dphi_);
  tmvaReader_->AddVariable("dlambda",&tmva_dlambda_);
  tmvaReader_->AddVariable("dqoverp",&tmva_dqoverp_);
  tmvaReader_->AddVariable("delta3d_r",&tmva_d3dr_);
  tmvaReader_->AddVariable("delta3d_z",&tmva_d3dz_);
  tmvaReader_->AddVariable("outer_nMissingInner",&tmva_outer_nMissingInner_);
  tmvaReader_->AddVariable("inner_nMissingOuter",&tmva_inner_nMissingOuter_);
  tmvaReader_->BookMVA("BDTG",mvaFilePath);
  */

  
}

DuplicateTrackMerger::~DuplicateTrackMerger()
{

  if(!useForestFromDB_) delete forest_;

}


#ifdef VI_STAT
  struct Stat {
    Stat() : maxCos(1.1), nCand(0),nLoop0(0) {}
    ~Stat() {
      std::cout << "Stats " << nCand << ' ' << nLoop0 << ' ' << maxCos << std::endl;
    }
    std::atomic<float> maxCos;
    std::atomic<int> nCand, nLoop0;
  };
  Stat stat;
#endif
  
template<typename T>
void update_maximum(std::atomic<T>& maximum_value, T const& value) noexcept
{
    T prev_value = maximum_value;
    while(prev_value < value &&
            !maximum_value.compare_exchange_weak(prev_value, value))
        ;
}

  template<typename T>
void update_minimum(std::atomic<T>& minimum_value, T const& value) noexcept
{
    T prev_value = minimum_value;
    while(prev_value > value &&
            !minimum_value.compare_exchange_weak(prev_value, value))
        ;
}

  
void DuplicateTrackMerger::produce(edm::Event& iEvent, const edm::EventSetup& iSetup)
{

  merger_.init(iSetup);

  if(!forest_){
    if(useForestFromDB_){
      edm::ESHandle<GBRForest> forestHandle;
      iSetup.get<GBRWrapperRcd>().get(forestLabel_,forestHandle);
      forest_ = forestHandle.product();
    }else{
      TFile gbrfile(dbFileName_.c_str());
      forest_ = dynamic_cast<const GBRForest*>(gbrfile.Get(forestLabel_.c_str()));
    }
  }

  //edm::Handle<edm::View<reco::Track> >handle;
  edm::Handle<reco::TrackCollection >handle;
  iEvent.getByToken(trackSource_,handle);
  auto const & tracks = *handle;
  
  edm::ESHandle<MagneticField> hmagfield;
  iSetup.get<IdealMagneticFieldRecord>().get(hmagfield);
  magfield_ = hmagfield.product();

  edm::ESHandle<TrackerTopology> httopo;
  iSetup.get<TrackerTopologyRcd>().get(httopo);
  ttopo_ = httopo.product();

  edm::ESHandle<TrackerGeometry> hgeom;
  iSetup.get<TrackerDigiGeometryRecord>().get(hgeom);
  geom_ = hgeom.product();

  edm::ESHandle<Propagator> hpropagator;
  iSetup.get<TrackingComponentsRecord>().get(propagatorName_, hpropagator);
  propagator_ = hpropagator.product();

  edm::ESHandle<Chi2MeasurementEstimatorBase> hestimator;
  iSetup.get<TrackingComponentsRecord>().get(chi2EstimatorName_, hestimator);
  chi2Estimator_ = hestimator.product();

  TSCPBuilderNoMaterial tscpBuilder;
  auto out_duplicateCandidates = std::make_unique<std::vector<TrackCandidate>>();

  auto out_candidateMap = std::make_unique<CandidateToDuplicate>();
  LogDebug("DuplicateTrackMerger") << "Number of tracks to be checked for merging: " << tracks.size();

#ifdef EDM_ML_DEBUG
  auto test = [&](const reco::Track *a, const reco::Track *b) {
    const auto ev = iEvent.id().event();
    const auto aOriAlgo = a->originalAlgo();
    const auto bOriAlgo = b->originalAlgo();
    const auto aSeed = a->seedRef().key();
    const auto bSeed = b->seedRef().key();
    return ((ev == 6903 && ((aOriAlgo == 23 && aSeed == 695 && bOriAlgo == 5 && bSeed == 652) ||
                            (aOriAlgo == 23 && aSeed == 400 && bOriAlgo == 7 && bSeed == 156) ||
                            (aOriAlgo == 4 && aSeed == 914 && bOriAlgo == 22 && bSeed == 503) ||
                            (aOriAlgo == 5 && aSeed == 809 && bOriAlgo ==  4 && bSeed == 1030) ||
                            (aOriAlgo == 23 && aSeed == 749 && bOriAlgo == 5 && bSeed == 659) ||
                            (aOriAlgo == 4 && aSeed == 1053 && bOriAlgo == 23 && bSeed == 1035) ||
                            (aOriAlgo == 4 && aSeed ==  810 && bOriAlgo == 5 && bSeed == 666) ||
                            (aOriAlgo == 4 && aSeed ==  974 && bOriAlgo == 5 && bSeed == 778))) ||
            (ev == 6904 && ((aOriAlgo == 23 && aSeed == 526 && bOriAlgo == 5 && bSeed == 307) ||
                            (aOriAlgo == 4 && aSeed == 559 && bOriAlgo == 22 && bSeed == 472))) ||
            (ev == 6902 && ((aOriAlgo == 4 && aSeed == 750 && bOriAlgo == 22 && bSeed == 340) ||
                            (aOriAlgo == 4 && aSeed == 906 && bOriAlgo == 5 && bSeed == 609) ||
                            (aOriAlgo == 4 && aSeed == 724 && bOriAlgo == 5 && bSeed == 528) ||
                            (aOriAlgo == 4 && aSeed == 943 && bOriAlgo == 23 && bSeed == 739) ||
                            (aOriAlgo == 8 && aSeed == 2 && bOriAlgo == 9 && bSeed == 2282) ||
                            (aOriAlgo == 23 && aSeed == 827 && bOriAlgo == 5 && bSeed == 656) ||
                            (aOriAlgo == 22 && aSeed == 667 && bOriAlgo == 7 && bSeed == 516))));
  };
#endif

  for(int i = 0; i < (int)tracks.size(); i++){
    const reco::Track *rt1 = &tracks[i];

    if(rt1->innerMomentum().perp2() < minpT2_)continue;
    // if(rt1->innerMomentum().R() < minP_)continue;
    for(int j = i+1; j < (int)tracks.size();j++){
      const reco::Track *rt2 = &tracks[j];

#ifdef EDM_ML_DEBUG
      debug_ = false;
      if(test(rt1, rt2) || test(rt2, rt1)) {
        debug_ = true;
        LogTrace("DuplicateTrackMerger") << "Track1 " << i << " originalAlgo " << rt1->originalAlgo() << " seed " << rt1->seedRef().key() << " pT " << std::sqrt(rt1->innerMomentum().perp2()) << " charge " << rt1->charge() << " outerPosition2 " << rt1->outerPosition().perp2() << "\n"
                                         << "Track2 " << j << " originalAlgo " << rt2->originalAlgo() << " seed " << rt2->seedRef().key() << " pT " << std::sqrt(rt2->innerMomentum().perp2()) << " charge " << rt2->charge() << " outerPosition2 " << rt2->outerPosition().perp2();
      }
#endif

      if(rt1->charge() != rt2->charge())continue;
      auto cosT = (*rt1).momentum().unit().Dot((*rt2).momentum().unit());
      IfLogTrace(debug_, "DuplicateTrackMerger") << " cosT " << cosT;
      if (cosT<0.) continue;
      if(rt2->innerMomentum().perp2() < minpT2_)continue;
      // if(rt2->innerMomentum().R() < minP_)continue;
      const reco::Track* t1,*t2;
      if(rt1->outerPosition().perp2() < rt2->outerPosition().perp2()){
    t1 = rt1;
    t2 = rt2;
      }else{
    t1 = rt2;
    t2 = rt1;
      }
      auto deltaR3d2 = (t1->outerPosition() - t2->innerPosition()).mag2();

      if(t1->outerPosition().perp2() > t2->innerPosition().perp2()) deltaR3d2 *= -1.0;
      IfLogTrace(debug_, "DuplicateTrackMerger") << " deltaR3d2 " << deltaR3d2 << " t1.outerPos2 " << t1->outerPosition().perp2() << " t2.innerPos2 " << t2->innerPosition().perp2();

      bool compatible = false;
      DuplicateTrackType duplType;
      if(deltaR3d2 >= minDeltaR3d2_) {
        compatible = checkForDisjointTracks(t1, t2, tscpBuilder);
        duplType = DuplicateTrackType::Disjoint;
      }
      else {
        compatible = checkForOverlappingTracks(t1, t2, cosT);
        duplType = DuplicateTrackType::Overlapping;
      }
      if(!compatible) continue;
      
      
      IfLogTrace(debug_, "DuplicateTrackMerger") << " marking as duplicates";
      out_duplicateCandidates->push_back(merger_.merge(*t1,*t2, duplType));
      out_candidateMap->emplace_back(i,j);

#ifdef VI_STAT
      ++stat.nCand;
      //    auto cosT = float((*t1).momentum().unit().Dot((*t2).momentum().unit()));
      if (cosT>0) update_minimum(stat.maxCos,float(cosT));
      else   ++stat.nLoop0;
#endif
      
    }
  }
  iEvent.put(std::move(out_duplicateCandidates),"candidates");
  iEvent.put(std::move(out_candidateMap),"candidateMap");

  
}


  bool DuplicateTrackMerger::checkForDisjointTracks(const reco::Track *t1, const reco::Track *t2, TSCPBuilderNoMaterial& tscpBuilder) const {
    IfLogTrace(debug_, "DuplicateTrackMerger") << " Checking for disjoint duplicates";

    FreeTrajectoryState fts1 = trajectoryStateTransform::outerFreeState(*t1, &*magfield_,false);
    FreeTrajectoryState fts2 = trajectoryStateTransform::innerFreeState(*t2, &*magfield_,false);
    GlobalPoint avgPoint((t1->outerPosition().x()+t2->innerPosition().x())*0.5,(t1->outerPosition().y()+t2->innerPosition().y())*0.5,(t1->outerPosition().z()+t2->innerPosition().z())*0.5);
    TrajectoryStateClosestToPoint TSCP1 = tscpBuilder(fts1, avgPoint);
    TrajectoryStateClosestToPoint TSCP2 = tscpBuilder(fts2, avgPoint);
    IfLogTrace(debug_, "DuplicateTrackMerger") << " TSCP1.isValid " << TSCP1.isValid() << " TSCP2.isValid " << TSCP2.isValid();
    if(!TSCP1.isValid()) return false;
    if(!TSCP2.isValid()) return false;

    const FreeTrajectoryState ftsn1 = TSCP1.theState();
    const FreeTrajectoryState ftsn2 = TSCP2.theState();
 
    IfLogTrace(debug_, "DuplicateTrackMerger") << " DCA2 " << (ftsn2.position()-ftsn1.position()).mag2();
    if ( (ftsn2.position()-ftsn1.position()).mag2() > maxDCA2_ ) return false;

    auto qoverp1 = ftsn1.signedInverseMomentum();
    auto qoverp2 = ftsn2.signedInverseMomentum();
    float tmva_dqoverp_ = qoverp1-qoverp2;
    IfLogTrace(debug_, "DuplicateTrackMerger") << " dqoverp " << tmva_dqoverp_;
    if ( std::abs(tmva_dqoverp_) > maxDQoP_ ) return false;


    //auto pp = [&](TrajectoryStateClosestToPoint const & ts) { std::cout << ' ' << ts.perigeeParameters().vector()[0] << '/'  << ts.perigeeError().transverseCurvatureError();};
    //if(qoverp1*qoverp2 <0) { std::cout << "charge different " << qoverp1 <<',' << qoverp2; pp(TSCP1); pp(TSCP2); std::cout << std::endl;}

    auto lambda1 =  M_PI/2 - ftsn1.momentum().theta();
    auto lambda2 =  M_PI/2 - ftsn2.momentum().theta();
    float tmva_dlambda_ = lambda1-lambda2;
    IfLogTrace(debug_, "DuplicateTrackMerger") << " dlambda " << tmva_dlambda_;
    if ( std::abs(tmva_dlambda_) > maxDLambda_ ) return false;

    auto phi1 = ftsn1.momentum().phi();
    auto phi2 = ftsn2.momentum().phi();
    float tmva_dphi_ = phi1-phi2;
    if(std::abs(tmva_dphi_) > float(M_PI)) tmva_dphi_ = 2.f*float(M_PI) - std::abs(tmva_dphi_);
    IfLogTrace(debug_, "DuplicateTrackMerger") << " dphi " << tmva_dphi_;
    if (std::abs(tmva_dphi_) > maxDPhi_ ) return false;

    auto dxy1 = (-ftsn1.position().x() * ftsn1.momentum().y() + ftsn1.position().y() * ftsn1.momentum().x())/TSCP1.pt();
    auto dxy2 = (-ftsn2.position().x() * ftsn2.momentum().y() + ftsn2.position().y() * ftsn2.momentum().x())/TSCP2.pt();
    float tmva_ddxy_ = dxy1-dxy2;
    IfLogTrace(debug_, "DuplicateTrackMerger") << " ddxy " << tmva_ddxy_;
    if ( std::abs(tmva_ddxy_) > maxDdxy_ ) return false;

    auto dsz1 = ftsn1.position().z() * TSCP1.pt() / TSCP1.momentum().mag()
      - (ftsn1.position().x() * ftsn1.momentum().y() + ftsn1.position().y() * ftsn1.momentum().x())/TSCP1.pt() * ftsn1.momentum().z()/ftsn1.momentum().mag();
    auto dsz2 = ftsn2.position().z() * TSCP2.pt() / TSCP2.momentum().mag()
      - (ftsn2.position().x() * ftsn2.momentum().y() + ftsn2.position().y() * ftsn2.momentum().x())/TSCP2.pt() * ftsn2.momentum().z()/ftsn2.momentum().mag();
    float tmva_ddsz_ = dsz1-dsz2;
    IfLogTrace(debug_, "DuplicateTrackMerger") << " ddsz " << tmva_ddsz_;
    if ( std::abs(tmva_ddsz_) > maxDdsz_ ) return false;

    float tmva_d3dr_ = avgPoint.perp();
    float tmva_d3dz_ = avgPoint.z();
    float tmva_outer_nMissingInner_ = t2->hitPattern().numberOfLostHits(reco::HitPattern::MISSING_INNER_HITS);
    float tmva_inner_nMissingOuter_ = t1->hitPattern().numberOfLostHits(reco::HitPattern::MISSING_OUTER_HITS);

    float gbrVals_[9];
    gbrVals_[0] = tmva_ddsz_;
    gbrVals_[1] = tmva_ddxy_;
    gbrVals_[2] = tmva_dphi_;
    gbrVals_[3] = tmva_dlambda_;
    gbrVals_[4] = tmva_dqoverp_;
    gbrVals_[5] = tmva_d3dr_;
    gbrVals_[6] = tmva_d3dz_;
    gbrVals_[7] = tmva_outer_nMissingInner_;
    gbrVals_[8] = tmva_inner_nMissingOuter_;

    auto mvaBDTG = forest_->GetClassifier(gbrVals_);
    IfLogTrace(debug_, "DuplicateTrackMerger") << " mvaBDTG " << mvaBDTG;
    if(mvaBDTG < minBDTG_) return false;

    //  std::cout << "to merge " << mvaBDTG << ' ' << std::copysign(std::sqrt(std::abs(deltaR3d2)),deltaR3d2) << ' ' << tmva_dphi_ << ' ' << TSCP1.pt() <<'/'<<TSCP2.pt() << std::endl;
    return true;
  }

  bool DuplicateTrackMerger::checkForOverlappingTracks(const reco::Track *t1, const reco::Track *t2, float cosT) const {
    // ensure t1 is the shorter track
    if(t2->numberOfValidHits() < t1->numberOfValidHits())
      std::swap(t1, t2);

    IfLogTrace(debug_, "DuplicateTrackMerger") << " Checking for overlapping duplicates, cosT " << cosT << " t1 hits " << t1->numberOfValidHits();
    if(cosT < overlapCheckMinCosT_) return false;
    if(t1->numberOfValidHits() > overlapCheckMaxHits_) return false;

    // find the hit on the longer track on layer of the first hit of the shorter track
    auto findHitOnT2 = [&](const TrackingRecHit *hit1) {
      const auto hitDet = hit1->geographicalId().det();
      const auto hitSubdet = hit1->geographicalId().subdetId();
      const auto hitLayer = ttopo_->layer(hit1->geographicalId());
      return std::find_if(t2->recHitsBegin(), t2->recHitsEnd(), [&](const TrackingRecHit *hit2) {
          const auto& detId = hit2->geographicalId();
          return (detId.det() == hitDet && detId.subdetId() == hitSubdet &&
                  ttopo_->layer(detId) == hitLayer);
        });
    };

    auto t1HitIter = t1->recHitsBegin();
    if(!(*t1HitIter)->isValid()) {
      IfLogTrace(debug_, "DuplicateTrackMerger") << " first t1 hit invalid";
      return false;
    }
    auto t2HitIter = findHitOnT2(*t1HitIter);
    if(t2HitIter == t2->recHitsEnd()) {
      // if first hit not found, try with second
      // if that fails, then reject
      ++t1HitIter;
      assert(t1HitIter != t1->recHitsEnd());

      if(!(*t1HitIter)->isValid()) {
        IfLogTrace(debug_, "DuplicateTrackMerger") << " second t1 hit invalid";
        return false;
      }
      t2HitIter = findHitOnT2(*t1HitIter);
      if(t2HitIter == t2->recHitsEnd()) return false;
    }
    IfLogTrace(debug_, "DuplicateTrackMerger") << " starting overlap check from t1 hit " << std::distance(t1->recHitsBegin(), t1HitIter)
                                               << " t2 hit " << std::distance(t2->recHitsBegin(), t2HitIter);

    auto prevSubdet = (*t1HitIter)->geographicalId().subdetId();
    const TrajectoryStateOnSurface tsosInner = trajectoryStateTransform::innerStateOnSurface(*t2, *geom_, magfield_);

    ++t1HitIter; ++t2HitIter;
    unsigned int missedLayers = 0;
    while(t1HitIter != t1->recHitsEnd() && t2HitIter != t2->recHitsEnd()) {
      // in case of invalid hits, reject immediately
      if((*t1HitIter)->getType() != TrackingRecHit::valid || trackerHitRTTI::isUndef(**t1HitIter) ||
         (*t2HitIter)->getType() != TrackingRecHit::valid || trackerHitRTTI::isUndef(**t2HitIter)) {
        IfLogTrace(debug_, "DuplicateTrackMerger") << " t1 hit " << std::distance(t1->recHitsBegin(), t1HitIter)
                                                   << " t2 hit " << std::distance(t2->recHitsBegin(), t2HitIter)
                                                   << " either is invalid, types t1 " << (*t1HitIter)->getType() << " t2 " << (*t2HitIter)->getType();
        return false;
      }

      const auto& t1DetId = (*t1HitIter)->geographicalId();
      const auto& t2DetId = (*t2HitIter)->geographicalId();

      const auto t1Det = t1DetId.det();
      const auto t2Det = t2DetId.det();
      if(t1Det != DetId::Tracker || t2Det != DetId::Tracker) {
          IfLogTrace(debug_, "DuplicateTrackMerger") << " t1 hit " << std::distance(t1->recHitsBegin(), t1HitIter)
                                                     << " t2 hit " << std::distance(t2->recHitsBegin(), t2HitIter)
                                                     << " either not from Tracker, dets t1 " << t1Det << " t2 " << t2Det;
        return false;
      }

      const auto t1Subdet = t1DetId.subdetId();
      const auto t1Layer = ttopo_->layer(t1DetId);

      // reject if hits have the same DetId but are different
      if(t1DetId == t2DetId) {
        if(!(*t1HitIter)->sharesInput(*t2HitIter, TrackingRecHit::all)) {
          IfLogTrace(debug_, "DuplicateTrackMerger") << " t1 hit " << std::distance(t1->recHitsBegin(), t1HitIter)
                                                     << " t2 hit " << std::distance(t2->recHitsBegin(), t2HitIter)
                                                     << " same DetId (" << t1DetId.rawId() << ") but do not share all input";
          return false;
        }
      }
      else {
        const auto t2Subdet = t2DetId.subdetId();
        const auto t2Layer = ttopo_->layer(t2DetId);

        // reject if hits are on different layers
        if(t1Subdet != t2Subdet || t1Layer != t2Layer) {
          bool recovered = false;
          // but try to recover first by checking if either one has skipped over a layer
          if(t1Subdet == prevSubdet && t2Subdet != prevSubdet) {
            // t1 has a layer t2 doesn't
            ++t1HitIter;
            recovered = true;
          }
          else if(t1Subdet != prevSubdet && t2Subdet == prevSubdet) {
            // t2 has a layer t1 doesn't
            ++t2HitIter;
            recovered = true;
          }
          else if(t1Subdet == t2Subdet) {
            prevSubdet = t1Subdet;
            // same subdet, so layer must be different
            if(t2Layer > t1Layer) {
              // t1 has a layer t2 doesn't
              ++t1HitIter;
              recovered = true;
            }
            else if(t1Layer > t2Layer) {
              // t2 has a layer t1 doesn't
              ++t2HitIter;
              recovered = true;
            }
          }
          if(recovered) {
            ++missedLayers;
            if(missedLayers > overlapCheckMaxMissingLayers_) {
              IfLogTrace(debug_, "DuplicateTrackMerger") << " max number of missed layers exceeded";
              return false;
            }
            continue;
          }

          IfLogTrace(debug_, "DuplicateTrackMerger") << " t1 hit " << std::distance(t1->recHitsBegin(), t1HitIter)
                                                     << " t2 hit " << std::distance(t2->recHitsBegin(), t2HitIter)
                                                     << " are on different layers (subdet, layer) t1 " << t1Subdet << "," << t1Layer
                                                     << " t2 " << t2Subdet << "," << t2Layer;
          return false;
        }
        // reject if same layer (but not same hit) in non-pixel detector
        else if(t1Subdet != PixelSubdetector::PixelBarrel && t1Subdet != PixelSubdetector::PixelEndcap) {
          IfLogTrace(debug_, "DuplicateTrackMerger") << " t1 hit " << std::distance(t1->recHitsBegin(), t1HitIter)
                                                     << " t2 hit " << std::distance(t2->recHitsBegin(), t2HitIter)
                                                     << " are on same layer, but in non-pixel detector (det " << t1Det << " subdet " << t1Subdet << " layer " << t1Layer << ")";
          return false;
        }
      }

      // Propagate longer track to the shorter track hit surface, check compatibility
      TrajectoryStateOnSurface tsosPropagated = propagator_->propagate(tsosInner, (*t1HitIter)->det()->surface());
      if(!tsosPropagated.isValid()) { // reject immediately if TSOS is not valid
        IfLogTrace(debug_, "DuplicateTrackMerger") << " t1 hit " << std::distance(t1->recHitsBegin(), t1HitIter)
                                                   << " t2 hit " << std::distance(t2->recHitsBegin(), t2HitIter)
                                                   << " TSOS not valid";
        return false;
      }
      auto passChi2Pair = chi2Estimator_->estimate(tsosPropagated, **t1HitIter);
      if(!passChi2Pair.first) {
        IfLogTrace(debug_, "DuplicateTrackMerger") << " t1 hit " << std::distance(t1->recHitsBegin(), t1HitIter)
                                                   << " t2 hit " << std::distance(t2->recHitsBegin(), t2HitIter)
                                                   << " hit chi2 compatibility failed with chi2 " << passChi2Pair.second;
        return false;
      }

      prevSubdet = t1Subdet;
      ++t1HitIter; ++t2HitIter;
    }
    if(t1HitIter != t1->recHitsEnd()) {
      IfLogTrace(debug_, "DuplicateTrackMerger") << " hits on t2 ended before hits on t1";
      return false;
    }

    IfLogTrace(debug_, "DuplicateTrackMerger") << " all hits on t2 are on layers whre t1 has also a hit";
    return true;
  }
}

#include "FWCore/PluginManager/interface/ModuleDef.h"
#include "FWCore/Framework/interface/MakerMacros.h"

DEFINE_FWK_MODULE(DuplicateTrackMerger);