PFElecTkProducer.cc

Go to the documentation of this file.

// -*- C++ -*-
//
// Package:    PFTracking
// Class:      PFElecTkProducer
// 
// Original Author:  Michele Pioppi
//         Created:  Tue Jan 23 15:26:39 CET 2007



// system include files
#include <memory>

// user include files
#include "RecoParticleFlow/PFTracking/interface/PFElecTkProducer.h"
#include "RecoParticleFlow/PFTracking/interface/PFTrackTransformer.h"
#include "RecoParticleFlow/PFTracking/interface/ConvBremPFTrackFinder.h"
#include "MagneticField/Engine/interface/MagneticField.h"
#include "MagneticField/Records/interface/IdealMagneticFieldRecord.h"
#include "FWCore/Framework/interface/ESHandle.h"
#include "TrackingTools/PatternTools/interface/Trajectory.h"
#include "DataFormats/TrajectorySeed/interface/TrajectorySeed.h"
#include "DataFormats/EgammaCandidates/interface/GsfElectronFwd.h"
#include "DataFormats/EgammaCandidates/interface/GsfElectron.h"
#include "Geometry/Records/interface/TrackerDigiGeometryRecord.h"
#include "Geometry/TrackerGeometryBuilder/interface/TrackerGeometry.h"
#include "MagneticField/Engine/interface/MagneticField.h"
#include "DataFormats/ParticleFlowReco/interface/PFClusterFwd.h"
#include "DataFormats/ParticleFlowReco/interface/PFCluster.h"
#include "TrackingTools/TransientTrack/interface/TransientTrackBuilder.h"
#include "TrackingTools/Records/interface/TransientTrackRecord.h"
#include "DataFormats/VertexReco/interface/Vertex.h"
#include "RecoParticleFlow/PFClusterTools/interface/LinkByRecHit.h"
#include "RecoParticleFlow/PFClusterTools/interface/ClusterClusterMapping.h"
#include "RecoParticleFlow/PFTracking/interface/PFTrackAlgoTools.h"

#include "TMath.h"
using namespace std;
using namespace edm;
using namespace reco;
PFElecTkProducer::PFElecTkProducer(const ParameterSet& iConfig, const convbremhelpers::HeavyObjectCache*):
  conf_(iConfig)
{
  

  gsfTrackLabel_ = consumes<reco::GsfTrackCollection>(iConfig.getParameter<InputTag>
                              ("GsfTrackModuleLabel"));

  pfTrackLabel_ = consumes<reco::PFRecTrackCollection>(iConfig.getParameter<InputTag>
                               ("PFRecTrackLabel"));

  primVtxLabel_ = consumes<reco::VertexCollection>(iConfig.getParameter<InputTag>
                           ("PrimaryVertexLabel"));

  pfEcalClusters_ = consumes<reco::PFClusterCollection>(iConfig.getParameter<InputTag>
                            ("PFEcalClusters"));

  pfNuclear_ = consumes<reco::PFDisplacedTrackerVertexCollection>(iConfig.getParameter<InputTag>
                                  ("PFNuclear"));
  
  pfConv_ = consumes<reco::PFConversionCollection>(iConfig.getParameter<InputTag>
                           ("PFConversions"));
  
  pfV0_ = consumes<reco::PFV0Collection>(iConfig.getParameter<InputTag>
                     ("PFV0"));

  useNuclear_ = iConfig.getParameter<bool>("useNuclear");
  useConversions_ = iConfig.getParameter<bool>("useConversions");
  useV0_ = iConfig.getParameter<bool>("useV0");
  debugGsfCleaning_ = iConfig.getParameter<bool>("debugGsfCleaning");

  produces<GsfPFRecTrackCollection>();
  produces<GsfPFRecTrackCollection>( "Secondary" ).setBranchAlias( "secondary" );


  trajinev_ = iConfig.getParameter<bool>("TrajInEvents");
  modemomentum_ = iConfig.getParameter<bool>("ModeMomentum");
  applySel_ = iConfig.getParameter<bool>("applyEGSelection");
  applyGsfClean_ = iConfig.getParameter<bool>("applyGsfTrackCleaning");
  applyAngularGsfClean_ = iConfig.getParameter<bool>("applyAlsoGsfAngularCleaning");
  detaCutGsfClean_ =  iConfig.getParameter<double>("maxDEtaGsfAngularCleaning");
  dphiCutGsfClean_ =  iConfig.getParameter<double>("maxDPhiBremTangGsfAngularCleaning");
  useFifthStepForTrackDriven_ = iConfig.getParameter<bool>("useFifthStepForTrackerDrivenGsf");
  useFifthStepForEcalDriven_ = iConfig.getParameter<bool>("useFifthStepForEcalDrivenGsf");
  maxPtConvReco_ = iConfig.getParameter<double>("MaxConvBremRecoPT");
  detaGsfSC_ = iConfig.getParameter<double>("MinDEtaGsfSC");
  dphiGsfSC_ = iConfig.getParameter<double>("MinDPhiGsfSC");
  SCEne_ = iConfig.getParameter<double>("MinSCEnergy");
  
  // set parameter for convBremFinder
  useConvBremFinder_ =     iConfig.getParameter<bool>("useConvBremFinder");
  
  mvaConvBremFinderIDBarrelLowPt_ = iConfig.getParameter<double>("pf_convBremFinderID_mvaCutBarrelLowPt");
  mvaConvBremFinderIDBarrelHighPt_ = iConfig.getParameter<double>("pf_convBremFinderID_mvaCutBarrelHighPt");
  mvaConvBremFinderIDEndcapsLowPt_ = iConfig.getParameter<double>("pf_convBremFinderID_mvaCutEndcapsLowPt");
  mvaConvBremFinderIDEndcapsHighPt_ = iConfig.getParameter<double>("pf_convBremFinderID_mvaCutEndcapsHighPt");
  
}


PFElecTkProducer::~PFElecTkProducer()
{ }


//
// member functions
//

// ------------ method called to produce the data  ------------
void
PFElecTkProducer::produce(Event& iEvent, const EventSetup& iSetup)
{


  //create the empty collections 
  auto gsfPFRecTrackCollection = std::make_unique<GsfPFRecTrackCollection>();
  
  auto gsfPFRecTrackCollectionSecondary = std::make_unique<GsfPFRecTrackCollection>();

  //read collections of tracks
  Handle<GsfTrackCollection> gsftrackscoll;
  iEvent.getByToken(gsfTrackLabel_,gsftrackscoll);

  //read collections of trajectories
  Handle<vector<Trajectory> > TrajectoryCollection;
 
  //read pfrectrack collection
  Handle<PFRecTrackCollection> thePfRecTrackCollection;
  iEvent.getByToken(pfTrackLabel_,thePfRecTrackCollection);
  const PFRecTrackCollection& PfRTkColl = *(thePfRecTrackCollection.product());

  // PFClusters
  Handle<PFClusterCollection> theECPfClustCollection;
  iEvent.getByToken(pfEcalClusters_,theECPfClustCollection);
  const PFClusterCollection& theEcalClusters = *(theECPfClustCollection.product());

  //Primary Vertexes
  Handle<reco::VertexCollection> thePrimaryVertexColl;
  iEvent.getByToken(primVtxLabel_,thePrimaryVertexColl);
  
  // Displaced Vertex
  Handle< reco::PFDisplacedTrackerVertexCollection > pfNuclears; 
  if( useNuclear_ ) 
    iEvent.getByToken(pfNuclear_, pfNuclears);
  
  // Conversions 
  Handle< reco::PFConversionCollection > pfConversions;
  if( useConversions_ ) 
    iEvent.getByToken(pfConv_,pfConversions);

  // V0
  Handle< reco::PFV0Collection > pfV0;
  if( useV0_ ) 
    iEvent.getByToken(pfV0_, pfV0);
    
      
  

  GsfTrackCollection gsftracks = *(gsftrackscoll.product());    
  vector<Trajectory> tjvec(0);
  if (trajinev_){
    iEvent.getByToken(gsfTrackLabel_,TrajectoryCollection); 
    
    tjvec= *(TrajectoryCollection.product());
  }
  
  
  vector<reco::GsfPFRecTrack> selGsfPFRecTracks;
  vector<reco::GsfPFRecTrack> primaryGsfPFRecTracks;
  std::map<unsigned int, std::vector<reco::GsfPFRecTrack> >  GsfPFMap;
  
  for( unsigned igsf = 0; igsf < gsftracks.size(); igsf++ ) {
    GsfTrackRef trackRef(gsftrackscoll, igsf);
    gsfInnerMomentumCache_.push_back(trackRef->pMode());
    TrajectoryStateOnSurface i_inTSOS = mtsTransform_.innerStateOnSurface((*trackRef));
    GlobalVector i_innMom;
    if(i_inTSOS.isValid()){
      multiTrajectoryStateMode::momentumFromModeCartesian(i_inTSOS,i_innMom);  
      gsfInnerMomentumCache_.back() = i_innMom.mag();
    }
  }

  for (unsigned int igsf=0; igsf<gsftracks.size();igsf++) {
    
    GsfTrackRef trackRef(gsftrackscoll, igsf);

    int kf_ind=FindPfRef(PfRTkColl,gsftracks[igsf],false);
    
    if (kf_ind>=0) {
      
      PFRecTrackRef kf_ref(thePfRecTrackCollection,
               kf_ind);

      // remove fifth step tracks
      if( useFifthStepForEcalDriven_ == false
      || useFifthStepForTrackDriven_ == false) {
    
    bool isFifthStepTrack = PFTrackAlgoTools::isFifthStep(kf_ref->trackRef()->algo());
    bool isEcalDriven = true;
    bool isTrackerDriven = true;
    
    if (trackRef->seedRef().get()==nullptr) {
      isEcalDriven = false;
      isTrackerDriven = false;
    }
    else {
      auto const& SeedFromRef= static_cast<ElectronSeed const&>(*(trackRef->extra()->seedRef()) );
      if(SeedFromRef.caloCluster().isNull())
        isEcalDriven = false;
      if(SeedFromRef.ctfTrack().isNull())
        isTrackerDriven = false;
    }
    //note: the same track could be both ecalDriven and trackerDriven
    if(isFifthStepTrack && 
       isEcalDriven &&
       isTrackerDriven == false &&
       useFifthStepForEcalDriven_ == false) {
      continue;
    }

    if(isFifthStepTrack && 
       isTrackerDriven  && 
       isEcalDriven == false && 
       useFifthStepForTrackDriven_ == false) {
      continue;
    }

    if(isFifthStepTrack && 
       isTrackerDriven && 
       isEcalDriven && 
       useFifthStepForTrackDriven_ == false &&
       useFifthStepForEcalDriven_ == false) {
      continue;
    }
      }
      
      pftrack_=GsfPFRecTrack( gsftracks[igsf].charge(), 
                  reco::PFRecTrack::GSF, 
                  igsf, trackRef,
                  kf_ref);
    } else  {
      PFRecTrackRef dummyRef;
      pftrack_=GsfPFRecTrack( gsftracks[igsf].charge(), 
                  reco::PFRecTrack::GSF, 
                  igsf, trackRef,
                  dummyRef);
    }
    
    
    bool validgsfbrem = false;
    if(trajinev_) {
      validgsfbrem = pfTransformer_->addPointsAndBrems(pftrack_, 
                               gsftracks[igsf], 
                               tjvec[igsf],
                               modemomentum_);
    } else {
      validgsfbrem = pfTransformer_->addPointsAndBrems(pftrack_, 
                               gsftracks[igsf], 
                               mtsTransform_);
    }
    
    bool passSel = true;
    if(applySel_) 
      passSel = applySelection(gsftracks[igsf]);      

    if(validgsfbrem && passSel) 
      selGsfPFRecTracks.push_back(pftrack_);
  }
  
  
  unsigned int count_primary = 0;
  if(!selGsfPFRecTracks.empty()) {
    for(unsigned int ipfgsf=0; ipfgsf<selGsfPFRecTracks.size();ipfgsf++) {
      
      vector<unsigned int> secondaries(0);
      secondaries.clear();
      bool keepGsf = true;

      if(applyGsfClean_) {
    keepGsf = resolveGsfTracks(selGsfPFRecTracks,ipfgsf,secondaries,theEcalClusters);
      }
 
      //is primary? 
      if(keepGsf == true) {

    // Find kf tracks from converted brem photons
    if(convBremFinder_->foundConvBremPFRecTrack(thePfRecTrackCollection,thePrimaryVertexColl,
                            pfNuclears,pfConversions,pfV0,
                                                    globalCache(),
                            useNuclear_,useConversions_,useV0_,
                            theEcalClusters,selGsfPFRecTracks[ipfgsf])) {
      const vector<PFRecTrackRef>& convBremPFRecTracks(convBremFinder_->getConvBremPFRecTracks());
      for(unsigned int ii = 0; ii<convBremPFRecTracks.size(); ii++) {
        selGsfPFRecTracks[ipfgsf].addConvBremPFRecTrackRef(convBremPFRecTracks[ii]);
      }
    }
    
    // save primaries gsf tracks
    //  gsfPFRecTrackCollection->push_back(selGsfPFRecTracks[ipfgsf]);
    primaryGsfPFRecTracks.push_back(selGsfPFRecTracks[ipfgsf]);
    
    
    // NOTE:: THE TRACKID IS USED TO LINK THE PRIMARY GSF TRACK. THIS NEEDS 
    // TO BE CHANGED AS SOON AS IT IS POSSIBLE TO CHANGE DATAFORMATS
    // A MODIFICATION HERE IMPLIES A MODIFICATION IN PFBLOCKALGO.CC/H
    unsigned int primGsfIndex = selGsfPFRecTracks[ipfgsf].trackId();
    vector<reco::GsfPFRecTrack> trueGsfPFRecTracks;
    if(!secondaries.empty()) {
      // loop on secondaries gsf tracks (from converted brems)
      for(unsigned int isecpfgsf=0; isecpfgsf<secondaries.size();isecpfgsf++) {
        
        PFRecTrackRef refsecKF =  selGsfPFRecTracks[(secondaries[isecpfgsf])].kfPFRecTrackRef();
        
        unsigned int secGsfIndex = selGsfPFRecTracks[(secondaries[isecpfgsf])].trackId();
        GsfTrackRef secGsfRef = selGsfPFRecTracks[(secondaries[isecpfgsf])].gsfTrackRef();

        if(refsecKF.isNonnull()) {
          // NOTE::IT SAVED THE TRACKID OF THE PRIMARY!!! THIS IS USED IN PFBLOCKALGO.CC/H
          secpftrack_= GsfPFRecTrack( gsftracks[secGsfIndex].charge(), 
                      reco::PFRecTrack::GSF, 
                      primGsfIndex, secGsfRef,
                      refsecKF);
        }
        else{
          PFRecTrackRef dummyRef;
          // NOTE::IT SAVED THE TRACKID OF THE PRIMARY!!! THIS IS USED IN PFBLOCKALGO.CC/H
          secpftrack_= GsfPFRecTrack( gsftracks[secGsfIndex].charge(), 
                      reco::PFRecTrack::GSF, 
                      primGsfIndex, secGsfRef,
                      dummyRef);
        }

        bool validgsfbrem = false;
        if(trajinev_) {
          validgsfbrem = pfTransformer_->addPointsAndBrems(secpftrack_,
                                   gsftracks[secGsfIndex], 
                                   tjvec[secGsfIndex],
                                   modemomentum_);
        } else {
          validgsfbrem = pfTransformer_->addPointsAndBrems(secpftrack_,
                                   gsftracks[secGsfIndex], 
                                   mtsTransform_);
        }         

        if(validgsfbrem) {
          gsfPFRecTrackCollectionSecondary->push_back(secpftrack_);
          trueGsfPFRecTracks.push_back(secpftrack_);
        }       
      }
    }
    GsfPFMap.insert(pair<unsigned int,std::vector<reco::GsfPFRecTrack> >(count_primary,trueGsfPFRecTracks));
    trueGsfPFRecTracks.clear();
    count_primary++;
      }
    }
  }
  
  
  const edm::OrphanHandle<GsfPFRecTrackCollection> gsfPfRefProd = 
    iEvent.put(std::move(gsfPFRecTrackCollectionSecondary),"Secondary");
  
  
  //now the secondary GsfPFRecTracks are in the event, the Ref can be created
  createGsfPFRecTrackRef(gsfPfRefProd,primaryGsfPFRecTracks,GsfPFMap);
  
  for(unsigned int iGSF = 0; iGSF<primaryGsfPFRecTracks.size();iGSF++){
    gsfPFRecTrackCollection->push_back(primaryGsfPFRecTracks[iGSF]);
  }
  iEvent.put(std::move(gsfPFRecTrackCollection));

  selGsfPFRecTracks.clear();
  GsfPFMap.clear();
  primaryGsfPFRecTracks.clear();

  std::vector<double>().swap(gsfInnerMomentumCache_);
}
  
// create the secondary GsfPFRecTracks Ref
void
PFElecTkProducer::createGsfPFRecTrackRef(const edm::OrphanHandle<reco::GsfPFRecTrackCollection>& gsfPfHandle,
                     std::vector<reco::GsfPFRecTrack>& gsfPFRecTrackPrimary,
                     const std::map<unsigned int, std::vector<reco::GsfPFRecTrack> >& MapPrimSec) {
  unsigned int cgsf=0;
  unsigned int csecgsf=0;
  for (std::map<unsigned int, std::vector<reco::GsfPFRecTrack> >::const_iterator igsf = MapPrimSec.begin();
       igsf != MapPrimSec.end(); igsf++,cgsf++) {
    vector<reco::GsfPFRecTrack> SecGsfPF = igsf->second;
    for (unsigned int iSecGsf=0; iSecGsf < SecGsfPF.size(); iSecGsf++) {
      edm::Ref<reco::GsfPFRecTrackCollection> refgprt(gsfPfHandle,csecgsf);
      gsfPFRecTrackPrimary[cgsf].addConvBremGsfPFRecTrackRef(refgprt);
      ++csecgsf;
    }
  }

  return;
}
// ------------- method for find the corresponding kf pfrectrack ---------------------
int
PFElecTkProducer::FindPfRef(const reco::PFRecTrackCollection  & PfRTkColl, 
                const reco::GsfTrack& gsftk,
                bool otherColl){


  if (gsftk.seedRef().get()==nullptr) return -1;
  auto const &  ElSeedFromRef=static_cast<ElectronSeed const&>( *(gsftk.extra()->seedRef()) );
  //CASE 1 ELECTRONSEED DOES NOT HAVE A REF TO THE CKFTRACK
  if (ElSeedFromRef.ctfTrack().isNull()){
    reco::PFRecTrackCollection::const_iterator pft=PfRTkColl.begin();
    reco::PFRecTrackCollection::const_iterator pftend=PfRTkColl.end();
    unsigned int i_pf=0;
    int ibest=-1;
    unsigned int ish_max=0;
    float dr_min=1000;
    //SEARCH THE PFRECTRACK THAT SHARES HITS WITH THE ELECTRON SEED
    // Here the cpu time can be improved. 
    for(;pft!=pftend;++pft){
      unsigned int ish=0;
      
      float dph= fabs(pft->trackRef()->phi()-gsftk.phi()); 
      if (dph>TMath::Pi()) dph-= TMath::TwoPi();
      float det=fabs(pft->trackRef()->eta()-gsftk.eta());
      float dr =sqrt(dph*dph+det*det);  
      
      for(auto const& hhit : pft->trackRef()->recHits()) {
        if (!hhit->isValid()) continue;
        TrajectorySeed::const_iterator hit = gsftk.seedRef()->recHits().first;
        TrajectorySeed::const_iterator hit_end = gsftk.seedRef()->recHits().second;
        for(;hit!=hit_end;++hit) {
          if (!(hit->isValid())) continue;
          if(hhit->sharesInput(&*(hit),TrackingRecHit::all)) ish++;
          // if((hit->geographicalId()==hhit->geographicalId())&&
          //     ((hhit->localPosition()-hit->localPosition()).mag()<0.01)) ish++;
        }
      }
      

      if ((ish>ish_max)||
      ((ish==ish_max)&&(dr<dr_min))){
    ish_max=ish;
    dr_min=dr;
    ibest=i_pf;
      }
      
   
    
      i_pf++;
    }
    if (ibest<0) return -1;
    
    if((ish_max==0) || (dr_min>0.05))return -1;
    if(otherColl && (ish_max==0)) return -1;
    return ibest;
  }
  else{
    //ELECTRON SEED HAS A REFERENCE
   
    reco::PFRecTrackCollection::const_iterator pft=PfRTkColl.begin();
    reco::PFRecTrackCollection::const_iterator pftend=PfRTkColl.end();
    unsigned int i_pf=0;
    
    for(;pft!=pftend;++pft){
      //REF COMPARISON
      if (pft->trackRef()==ElSeedFromRef.ctfTrack()){
    return i_pf;
      }
      i_pf++;
    }
  }
  return -1;
}

// -- method to apply gsf electron selection to EcalDriven seeds
bool 
PFElecTkProducer::applySelection(const reco::GsfTrack& gsftk) {
  if (gsftk.seedRef().get()==nullptr) return false;
  auto const& ElSeedFromRef=static_cast<ElectronSeed const&>( *(gsftk.extra()->seedRef()) );

  bool passCut = false;
  if (ElSeedFromRef.ctfTrack().isNull()){
    if(ElSeedFromRef.caloCluster().isNull()) return passCut;
    auto const* scRef = static_cast<SuperCluster const*>(ElSeedFromRef.caloCluster().get());
    //do this just to know if exist a SC? 
    if(scRef) {
      float caloEne = scRef->energy();
      float feta = fabs(scRef->eta()-gsftk.etaMode());
      float fphi = fabs(scRef->phi()-gsftk.phiMode());
      if (fphi>TMath::Pi()) fphi-= TMath::TwoPi();
      if(caloEne > SCEne_ && feta < detaGsfSC_ && fabs(fphi) < dphiGsfSC_)
    passCut = true;
    }
  }
  else {
    // get all the gsf found by tracker driven
    passCut = true;
  }
  return passCut;
}
bool 
PFElecTkProducer::resolveGsfTracks(const vector<reco::GsfPFRecTrack>  & GsfPFVec, 
                   unsigned int ngsf, 
                   vector<unsigned int> &secondaries,
                   const reco::PFClusterCollection & theEClus) {
  bool debugCleaning = debugGsfCleaning_;
  bool n_keepGsf = true;

  const reco::GsfTrackRef& nGsfTrack = GsfPFVec[ngsf].gsfTrackRef();
  
  if (nGsfTrack->seedRef().get()==nullptr) return false;    
  auto const& nElSeedFromRef=static_cast<ElectronSeed const&>( *(nGsfTrack->extra()->seedRef()) );
  
  /* // now gotten from cache below
  TrajectoryStateOnSurface inTSOS = mtsTransform_.innerStateOnSurface((*nGsfTrack));
  GlobalVector ninnMom;
  float nPin =  nGsfTrack->pMode();
  if(inTSOS.isValid()){
    multiTrajectoryStateMode::momentumFromModeCartesian(inTSOS,ninnMom);
    nPin = ninnMom.mag();
  }
  */
  float nPin = gsfInnerMomentumCache_[nGsfTrack.key()];

  float neta = nGsfTrack->innerMomentum().eta();
  float nphi = nGsfTrack->innerMomentum().phi();
  


  
  if(debugCleaning)
    cout << " PFElecTkProducer:: considering track " << nGsfTrack->pt() 
     << " eta,phi " <<  nGsfTrack->eta() << ", " <<  nGsfTrack->phi()  << endl;
  
  
  for (unsigned int igsf=0; igsf< GsfPFVec.size();igsf++) {
    if(igsf != ngsf ) {
      reco::GsfTrackRef iGsfTrack = GsfPFVec[igsf].gsfTrackRef();

      if(debugCleaning)
    cout << " PFElecTkProducer:: and  comparing with track " << iGsfTrack->pt() 
         << " eta,phi " <<  iGsfTrack->eta() << ", " <<  iGsfTrack->phi()  << endl;
      
      float ieta = iGsfTrack->innerMomentum().eta();
      float iphi = iGsfTrack->innerMomentum().phi();
      float feta = fabs(neta - ieta);
      float fphi = fabs(nphi - iphi);
      if (fphi>TMath::Pi()) fphi-= TMath::TwoPi();     
  

      // apply a superloose preselection only to avoid un-useful cpu time: hard-coded for this reason
      if(feta < 0.5 && fabs(fphi) < 1.0) {
    if(debugCleaning)
      cout << " Entering angular superloose preselection " << endl;
    
        /* //now taken from cache below
    TrajectoryStateOnSurface i_inTSOS = mtsTransform_.innerStateOnSurface((*iGsfTrack));
    GlobalVector i_innMom;
    float iPin = iGsfTrack->pMode();
    if(i_inTSOS.isValid()){
      multiTrajectoryStateMode::momentumFromModeCartesian(i_inTSOS,i_innMom);  
      iPin = i_innMom.mag();
    }
        */
        float iPin = gsfInnerMomentumCache_[iGsfTrack.key()];

    if (iGsfTrack->seedRef().get()==nullptr) continue;   
    auto const& iElSeedFromRef=static_cast<ElectronSeed const&>( *(iGsfTrack->extra()->seedRef()) );

    float SCEnergy = -1.;
    // Check if two tracks match the same SC     
    bool areBothGsfEcalDriven = false;;
    bool isSameSC = isSameEgSC(nElSeedFromRef,iElSeedFromRef,areBothGsfEcalDriven,SCEnergy);
    
    // CASE1 both GsfTracks ecalDriven and match the same SC
    if(areBothGsfEcalDriven ) {
      if(isSameSC) {
        float nEP = SCEnergy/nPin;
        float iEP =  SCEnergy/iPin;
        if(debugCleaning)
          cout << " Entering SAME supercluster case " 
           << " nEP " << nEP 
           << " iEP " << iEP << endl;
        
        
        
        // if same SC take the closest or if same
        // distance the best E/p
        
        // Innermost using LostHits technology
        bool isSameLayer = false;
        bool iGsfInnermostWithLostHits = 
          isInnerMostWithLostHits(nGsfTrack,iGsfTrack,isSameLayer);
        
        
        if(debugCleaning)
          cout << " iGsf is InnerMostWithLostHits " << iGsfInnermostWithLostHits
           << " isSameLayer " << isSameLayer  << endl;
        
        if (iGsfInnermostWithLostHits) {
          n_keepGsf = false;
          return n_keepGsf;
        }
        else if(isSameLayer){
          if(fabs(iEP-1) < fabs(nEP-1)) {
        n_keepGsf = false;
        return n_keepGsf;
          }
          else{
        secondaries.push_back(igsf);
          }
        }
        else {
          // save secondaries gsf track (put selection)
        secondaries.push_back(igsf);
        }
      } // end same SC case
    }
    else {
      // enter in the condition where at least one track is trackerDriven
      float minBremDphi =  minTangDist(GsfPFVec[ngsf],GsfPFVec[igsf]);
      float nETot = 0.;
      float iETot = 0.;
      bool isBothGsfTrackerDriven = false;
      bool nEcalDriven = false;
      bool iEcalDriven = false;
      bool isSameScEgPf = isSharingEcalEnergyWithEgSC(GsfPFVec[ngsf],
                              GsfPFVec[igsf],
                              nElSeedFromRef,
                              iElSeedFromRef,
                              theEClus,
                              isBothGsfTrackerDriven,
                              nEcalDriven,
                              iEcalDriven,
                              nETot,
                              iETot);

      // check if the first hit of iGsfTrack < nGsfTrack          
      bool isSameLayer = false;
      bool iGsfInnermostWithLostHits = 
        isInnerMostWithLostHits(nGsfTrack,iGsfTrack,isSameLayer);

      if(isSameScEgPf) {
        // CASE 2 : One Gsf has reference to a SC and the other one not or both not
            
        if(debugCleaning) {
          cout << " Sharing ECAL energy passed " 
           << " nEtot " << nETot 
           << " iEtot " << iETot << endl;
          if(isBothGsfTrackerDriven) 
        cout << " Both Track are trackerDriven " << endl;
        }

        // Innermost using LostHits technology
        if (iGsfInnermostWithLostHits) {
          n_keepGsf = false;
          return n_keepGsf;
        }
        else if(isSameLayer){
          // Thirt Case:  One Gsf has reference to a SC and the other one not or both not
          // gsf tracks starts from the same layer
          // check number of sharing modules (at least 50%)
          // check number of sharing hits (at least 2)
          // check charge flip inner/outer

         
        // they share energy
          if(isBothGsfTrackerDriven == false) {
        // if at least one Gsf track is EcalDriven choose that one.
        if(iEcalDriven) {
          n_keepGsf = false;
          return n_keepGsf;
        }
        else {
          secondaries.push_back(igsf);
        }
          }
          else {
        // if both tracks are tracker driven choose that one with the best E/p
        // with ETot = max(En,Ei)

        float ETot = -1;
        if(nETot != iETot) {
          if(nETot > iETot)
            ETot = nETot;
          else 
            ETot = iETot;
        }
        else {
          ETot = nETot;
        }
        float nEP = ETot/nPin;
        float iEP = ETot/iPin;
        
        
        if(debugCleaning) 
          cout << " nETot " << nETot
               << " iETot " << iETot 
               << " ETot " << ETot << endl 
               << " nPin " << nPin
               << " iPin " << iPin 
               << " nEP " << nEP 
               << " iEP " << iEP << endl;
        
    
        if(fabs(iEP-1) < fabs(nEP-1)) {
          n_keepGsf = false;
          return n_keepGsf;
        }
        else{
          secondaries.push_back(igsf);
        }
          }
        }
        else {
          secondaries.push_back(igsf);
        }
      }
      else if(feta < detaCutGsfClean_ && minBremDphi < dphiCutGsfClean_) {
        // very close tracks
        bool secPushedBack = false;
        if(nEcalDriven == false && nETot == 0.) {
          n_keepGsf = false;
          return n_keepGsf;
        }
        else if(iEcalDriven == false && iETot == 0.) {
          secondaries.push_back(igsf);
          secPushedBack = true;
        }
        if(debugCleaning)
          cout << " Close Tracks " 
           << " feta " << feta << " fabs(fphi) " << fabs(fphi) 
           << " minBremDphi " <<  minBremDphi 
           << " nETot " << nETot 
           << " iETot " << iETot 
           << " nLostHits " <<  nGsfTrack->hitPattern().numberOfLostHits(HitPattern::MISSING_INNER_HITS) 
           << " iLostHits " << iGsfTrack->hitPattern().numberOfLostHits(HitPattern::MISSING_INNER_HITS) << endl;
        
        // apply selection only if one track has lost hits
        if(applyAngularGsfClean_) {
          if (iGsfInnermostWithLostHits) {
        n_keepGsf = false;
        return n_keepGsf;
          }
          else if(isSameLayer == false) {
        if(secPushedBack == false) 
          secondaries.push_back(igsf);
          }
        }
      }
      else if(feta < 0.1 && minBremDphi < 0.2){
        // failed all the conditions, discard only tracker driven tracks
        // with no PFClusters linked. 
        if(debugCleaning)
          cout << " Close Tracks and failed all the conditions " 
           << " feta " << feta << " fabs(fphi) " << fabs(fphi) 
           << " minBremDphi " <<  minBremDphi 
           << " nETot " << nETot 
           << " iETot " << iETot 
           << " nLostHits " <<  nGsfTrack->hitPattern().numberOfLostHits(HitPattern::MISSING_INNER_HITS) 
           << " iLostHits " << iGsfTrack->hitPattern().numberOfLostHits(HitPattern::MISSING_INNER_HITS) << endl;
        
        if(nEcalDriven == false && nETot == 0.) {
          n_keepGsf = false;
          return n_keepGsf;
        }
        // Here I do not push back the secondary because considered fakes...
      }
    }
      }
    }
  }
  
  return n_keepGsf;
}
float 
PFElecTkProducer::minTangDist(const reco::GsfPFRecTrack& primGsf,
                  const reco::GsfPFRecTrack& secGsf) {

  float minDphi = 1000.; 


  std::vector<reco::PFBrem> primPFBrem = primGsf.PFRecBrem();
  std::vector<reco::PFBrem> secPFBrem = secGsf.PFRecBrem();


  unsigned int cbrem = 0;
  for (unsigned isbrem = 0; isbrem < secPFBrem.size(); isbrem++) {
    if(secPFBrem[isbrem].indTrajPoint() == 99) continue;
    const reco::PFTrajectoryPoint& atSecECAL 
      = secPFBrem[isbrem].extrapolatedPoint( reco::PFTrajectoryPoint::ECALEntrance );
    if( ! atSecECAL.isValid() ) continue;
    float secPhi  = atSecECAL.positionREP().Phi();

    unsigned int sbrem = 0;
    for(unsigned ipbrem = 0; ipbrem < primPFBrem.size(); ipbrem++) {
      if(primPFBrem[ipbrem].indTrajPoint() == 99) continue;
      const reco::PFTrajectoryPoint& atPrimECAL 
    = primPFBrem[ipbrem].extrapolatedPoint( reco::PFTrajectoryPoint::ECALEntrance );
      if( ! atPrimECAL.isValid() ) continue;
      sbrem++;
      if(sbrem <= 3) {
    float primPhi = atPrimECAL.positionREP().Phi();
    
    float dphi = fabs(primPhi - secPhi);
    if (dphi>TMath::Pi()) dphi-= TMath::TwoPi();     
    if(fabs(dphi) < minDphi) {     
      minDphi = fabs(dphi);
    }
      }
    }
    
    
    cbrem++;
    if(cbrem == 3) 
      break;
  }
  return minDphi;
}
bool 
PFElecTkProducer::isSameEgSC(const reco::ElectronSeed& nSeed,
                 const reco::ElectronSeed& iSeed,
                 bool& bothGsfEcalDriven,
                 float& SCEnergy) {
  
  bool isSameSC = false;

  if(nSeed.caloCluster().isNonnull() && iSeed.caloCluster().isNonnull()) {
    auto const* nscRef = static_cast<SuperCluster const*>(nSeed.caloCluster().get());
    auto const* iscRef = static_cast<SuperCluster const*>(iSeed.caloCluster().get());

    if(nscRef && iscRef) {
      bothGsfEcalDriven = true;
      if(nscRef == iscRef) {
    isSameSC = true;
    // retrieve the supercluster energy
    SCEnergy = nscRef->energy();
      }
    }
  }
  return isSameSC;
}
bool 
PFElecTkProducer::isSharingEcalEnergyWithEgSC(const reco::GsfPFRecTrack& nGsfPFRecTrack,
                          const reco::GsfPFRecTrack& iGsfPFRecTrack,
                          const reco::ElectronSeed& nSeed,
                          const reco::ElectronSeed& iSeed,
                          const reco::PFClusterCollection& theEClus,
                          bool& bothGsfTrackerDriven,
                          bool& nEcalDriven,
                          bool& iEcalDriven,
                          float& nEnergy,
                          float& iEnergy) {
  
  bool isSharingEnergy = false;

  //which is EcalDriven?
  bool oneEcalDriven = true;
  SuperCluster const* scRef = nullptr;
  GsfPFRecTrack gsfPfTrack;

  if(nSeed.caloCluster().isNonnull()) {
    scRef = static_cast<SuperCluster const*>(nSeed.caloCluster().get());
    assert(scRef);
    nEnergy = scRef->energy();
    nEcalDriven = true;
    gsfPfTrack = iGsfPFRecTrack;
  }
  else if(iSeed.caloCluster().isNonnull()){
    scRef = static_cast<SuperCluster const*>(iSeed.caloCluster().get());
    assert(scRef);
    iEnergy = scRef->energy();
    iEcalDriven = true;
    gsfPfTrack = nGsfPFRecTrack;
  }
  else{
     oneEcalDriven = false;
  }
  
  if(oneEcalDriven) {
    //run a basic reconstruction for the particle flow

    vector<PFCluster> vecPFClusters;
    vecPFClusters.clear();

    for (PFClusterCollection::const_iterator clus = theEClus.begin();
     clus != theEClus.end();
     clus++ ) {
      PFCluster clust = *clus;
      clust.calculatePositionREP();

      float deta = fabs(scRef->position().eta() - clust.position().eta());
      float dphi = fabs(scRef->position().phi() - clust.position().phi());
      if (dphi>TMath::Pi()) dphi-= TMath::TwoPi();
 
      // Angle preselection between the supercluster and pfclusters
      // this is needed just to save some cpu-time for this is hard-coded     
      if(deta < 0.5 && fabs(dphi) < 1.0) {
    double distGsf = gsfPfTrack.extrapolatedPoint( reco::PFTrajectoryPoint::ECALShowerMax ).isValid() ?
      LinkByRecHit::testTrackAndClusterByRecHit(gsfPfTrack , clust ) : -1.;
    // check if it touch the GsfTrack
    if(distGsf > 0.) {
      if(nEcalDriven) 
        iEnergy += clust.energy();
      else
        nEnergy += clust.energy();
      vecPFClusters.push_back(clust);
    }
    // check if it touch the Brem-tangents
    else {
      vector<PFBrem> primPFBrem = gsfPfTrack.PFRecBrem();
      for(unsigned ipbrem = 0; ipbrem < primPFBrem.size(); ipbrem++) {
        if(primPFBrem[ipbrem].indTrajPoint() == 99) continue;
        const reco::PFRecTrack& pfBremTrack = primPFBrem[ipbrem];
        double dist = pfBremTrack.extrapolatedPoint( reco::PFTrajectoryPoint::ECALShowerMax ).isValid() ?
          LinkByRecHit::testTrackAndClusterByRecHit(pfBremTrack , clust, true ) : -1.;
        if(dist > 0.) {
          if(nEcalDriven) 
        iEnergy += clust.energy();
          else
        nEnergy += clust.energy();
          vecPFClusters.push_back(clust);
        }
      } 
    }  
      } // END if angle preselection
    } // PFClusters Loop
    if(!vecPFClusters.empty() ) {
      for(unsigned int pf = 0; pf < vecPFClusters.size(); pf++) {
    bool isCommon = ClusterClusterMapping::overlap(vecPFClusters[pf],*scRef);
    if(isCommon) {
      isSharingEnergy = true;
    }
    break;
      }
    }
  }
  else {
    // both tracks are trackerDriven, try ECAL energy matching also in this case.

    bothGsfTrackerDriven = true;
    vector<PFCluster> nPFCluster;
    vector<PFCluster> iPFCluster;

    nPFCluster.clear();
    iPFCluster.clear();

    for (PFClusterCollection::const_iterator clus = theEClus.begin();
     clus != theEClus.end();
     clus++ ) {
      PFCluster clust = *clus;
      clust.calculatePositionREP();
      
      float ndeta = fabs(nGsfPFRecTrack.gsfTrackRef()->eta() - clust.position().eta());
      float ndphi = fabs(nGsfPFRecTrack.gsfTrackRef()->phi() - clust.position().phi());
      if (ndphi>TMath::Pi()) ndphi-= TMath::TwoPi();
      // Apply loose preselection with the track
      // just to save cpu time, for this hard-coded
      if(ndeta < 0.5 && fabs(ndphi) < 1.0) {
    
    double distGsf = nGsfPFRecTrack.extrapolatedPoint( reco::PFTrajectoryPoint::ECALShowerMax ).isValid() ?
      LinkByRecHit::testTrackAndClusterByRecHit(nGsfPFRecTrack , clust ) : -1.;
    if(distGsf > 0.) {
      nPFCluster.push_back(clust);
      nEnergy += clust.energy();
    }
    else {
      const vector<PFBrem>& primPFBrem = nGsfPFRecTrack.PFRecBrem();
      for(unsigned ipbrem = 0; ipbrem < primPFBrem.size(); ipbrem++) {
        if(primPFBrem[ipbrem].indTrajPoint() == 99) continue;
        const reco::PFRecTrack& pfBremTrack = primPFBrem[ipbrem];
        double dist = pfBremTrack.extrapolatedPoint( reco::PFTrajectoryPoint::ECALShowerMax ).isValid() ?
        LinkByRecHit::testTrackAndClusterByRecHit(pfBremTrack , clust, true ) : -1.;
        if(dist > 0.) {
          nPFCluster.push_back(clust);
          nEnergy += clust.energy();
          break;
        }
      }
    }
      }

      float ideta = fabs(iGsfPFRecTrack.gsfTrackRef()->eta() - clust.position().eta());
      float idphi = fabs(iGsfPFRecTrack.gsfTrackRef()->phi() - clust.position().phi());
      if (idphi>TMath::Pi()) idphi-= TMath::TwoPi();
      // Apply loose preselection with the track
      // just to save cpu time, for this hard-coded
      if(ideta < 0.5 && fabs(idphi) < 1.0) {
    double dist = iGsfPFRecTrack.extrapolatedPoint( reco::PFTrajectoryPoint::ECALShowerMax ).isValid() ?
      LinkByRecHit::testTrackAndClusterByRecHit(iGsfPFRecTrack , clust ) : -1.;
    if(dist > 0.) {
      iPFCluster.push_back(clust);
      iEnergy += clust.energy();
    }
    else {
      vector<PFBrem> primPFBrem = iGsfPFRecTrack.PFRecBrem();
      for(unsigned ipbrem = 0; ipbrem < primPFBrem.size(); ipbrem++) {
        if(primPFBrem[ipbrem].indTrajPoint() == 99) continue;
        const reco::PFRecTrack& pfBremTrack = primPFBrem[ipbrem];
        double dist = LinkByRecHit::testTrackAndClusterByRecHit(pfBremTrack , clust, true );
        if(dist > 0.) {
          iPFCluster.push_back(clust);
          iEnergy += clust.energy();
          break;
        }
      }
    }
      }
    }


    if(!nPFCluster.empty() && !iPFCluster.empty()) {
      for(unsigned int npf = 0; npf < nPFCluster.size(); npf++) {
    for(unsigned int ipf = 0; ipf < iPFCluster.size(); ipf++) {
      bool isCommon = ClusterClusterMapping::overlap(nPFCluster[npf],iPFCluster[ipf]);
      if(isCommon) {
        isSharingEnergy = true;
        break;
      }
    }
    if(isSharingEnergy)
      break;
      }
    }
  }

  return isSharingEnergy;
}
bool PFElecTkProducer::isInnerMost(const reco::GsfTrackRef& nGsfTrack,
                   const reco::GsfTrackRef& iGsfTrack,
                   bool& sameLayer) {
  
  // copied by the class RecoEgamma/EgammaElectronAlgos/src/EgAmbiguityTools.cc
  // obsolete but the code is kept: now using lost hits method

  const reco::HitPattern &gsfHitPattern1 = nGsfTrack->hitPattern();
  const reco::HitPattern &gsfHitPattern2 = iGsfTrack->hitPattern();
  
  // retrieve first valid hit
  int gsfHitCounter1 = 0 ;
  for(auto const& hit : nGsfTrack->recHits()) {
    if (hit->isValid()) break ;
    gsfHitCounter1++;
  }
  
  int gsfHitCounter2 = 0 ;
  for(auto const& hit : iGsfTrack->recHits()) {
    if (hit->isValid()) break ;
    gsfHitCounter2++;
  }
  
  uint32_t gsfHit1 = gsfHitPattern1.getHitPattern(HitPattern::TRACK_HITS, gsfHitCounter1) ;
  uint32_t gsfHit2 = gsfHitPattern2.getHitPattern(HitPattern::TRACK_HITS, gsfHitCounter2) ;
  
  
  if (gsfHitPattern1.getSubStructure(gsfHit1)!=gsfHitPattern2.getSubStructure(gsfHit2))
   { 
     return (gsfHitPattern2.getSubStructure(gsfHit2)<gsfHitPattern1.getSubStructure(gsfHit1)); 
   }
  else if (gsfHitPattern1.getLayer(gsfHit1)!=gsfHitPattern2.getLayer(gsfHit2))
    { 
      return (gsfHitPattern2.getLayer(gsfHit2)<gsfHitPattern1.getLayer(gsfHit1)); 
    }
  else
   { 
     sameLayer = true;
     return  false; 
   }
}
bool PFElecTkProducer::isInnerMostWithLostHits(const reco::GsfTrackRef& nGsfTrack,
                           const reco::GsfTrackRef& iGsfTrack,
                           bool& sameLayer) {
  
  // define closest using the lost hits on the expectedhitsineer
  unsigned int nLostHits = nGsfTrack->hitPattern().numberOfLostHits(HitPattern::MISSING_INNER_HITS);
  unsigned int iLostHits = iGsfTrack->hitPattern().numberOfLostHits(HitPattern::MISSING_INNER_HITS);
  
  if (nLostHits!=iLostHits) {
    return (nLostHits > iLostHits);
  } 
  else {
    sameLayer = true;
    return  false; 
  }
}



// ------------ method called once each job just before starting event loop  ------------
void 
PFElecTkProducer::beginRun(const edm::Run& run,
               const EventSetup& iSetup)
{
  ESHandle<MagneticField> magneticField;
  iSetup.get<IdealMagneticFieldRecord>().get(magneticField);

  ESHandle<TrackerGeometry> tracker;
  iSetup.get<TrackerDigiGeometryRecord>().get(tracker);

  mtsTransform_ = MultiTrajectoryStateTransform(tracker.product(),magneticField.product());  

  pfTransformer_.reset( new PFTrackTransformer(math::XYZVector(magneticField->inTesla(GlobalPoint(0,0,0)))) );  
  
  edm::ESHandle<TransientTrackBuilder> builder;
  iSetup.get<TransientTrackRecord>().get("TransientTrackBuilder", builder);
  TransientTrackBuilder thebuilder = *(builder.product());
  
  convBremFinder_.reset( new ConvBremPFTrackFinder(thebuilder,
                                                   mvaConvBremFinderIDBarrelLowPt_,
                                                   mvaConvBremFinderIDBarrelHighPt_,
                                                   mvaConvBremFinderIDEndcapsLowPt_, 
                                                   mvaConvBremFinderIDEndcapsHighPt_) );
 
}

// ------------ method called once each job just after ending the event loop  ------------
void 
PFElecTkProducer::endRun(const edm::Run& run,
             const EventSetup& iSetup) {
  pfTransformer_.reset();
  convBremFinder_.reset();
}

//define this as a plug-in