PFTracking/plugins/PFElecTkProducer.cc

0001 // -*- C++ -*-
0002 //
0003 // Package:    PFTracking
0004 // Class:      PFElecTkProducer
0005 //
0006 // Original Author:  Michele Pioppi
0007 //         Created:  Tue Jan 23 15:26:39 CET 2007
0008
0009 /// \brief Abstract
0010 /*!
0011 \author Michele Pioppi, Daniele Benedetti
0012 \date January 2007
0013
0014  PFElecTkProducer reads the merged GsfTracks collection
0015  built with the TrackerDriven and EcalDriven seeds
0016  and transform them in PFGsfRecTracks.
0017 */
0018
0019 #include <memory>
0020
0021 #include "CommonTools/Utils/interface/KinematicTables.h"
0022 #include "CommonTools/Utils/interface/LazyConstructed.h"
0023 #include "DataFormats/EgammaCandidates/interface/GsfElectron.h"
0024 #include "DataFormats/EgammaReco/interface/ElectronSeed.h"
0025 #include "DataFormats/GsfTrackReco/interface/GsfTrack.h"
0026 #include "DataFormats/Math/interface/deltaR.h"
0027 #include "DataFormats/ParticleFlowReco/interface/GsfPFRecTrack.h"
0028 #include "DataFormats/ParticleFlowReco/interface/PFCluster.h"
0029 #include "DataFormats/ParticleFlowReco/interface/PFConversion.h"
0030 #include "DataFormats/ParticleFlowReco/interface/PFDisplacedTrackerVertex.h"
0031 #include "DataFormats/ParticleFlowReco/interface/PFDisplacedVertex.h"
0032 #include "DataFormats/ParticleFlowReco/interface/PFRecTrack.h"
0033 #include "DataFormats/ParticleFlowReco/interface/PFV0.h"
0034 #include "DataFormats/TrajectorySeed/interface/TrajectorySeed.h"
0035 #include "DataFormats/VertexReco/interface/Vertex.h"
0036 #include "FWCore/Framework/interface/ESHandle.h"
0037 #include "FWCore/Framework/interface/Event.h"
0038 #include "FWCore/Framework/interface/stream/EDProducer.h"
0039 #include "FWCore/ParameterSet/interface/ParameterSet.h"
0040 #include "Geometry/Records/interface/TrackerDigiGeometryRecord.h"
0041 #include "Geometry/TrackerGeometryBuilder/interface/TrackerGeometry.h"
0042 #include "MagneticField/Engine/interface/MagneticField.h"
0043 #include "MagneticField/Records/interface/IdealMagneticFieldRecord.h"
0044 #include "RecoParticleFlow/PFClusterTools/interface/ClusterClusterMapping.h"
0045 #include "RecoParticleFlow/PFClusterTools/interface/LinkByRecHit.h"
0046 #include "RecoParticleFlow/PFTracking/interface/ConvBremHeavyObjectCache.h"
0047 #include "RecoParticleFlow/PFTracking/interface/ConvBremPFTrackFinder.h"
0048 #include "RecoParticleFlow/PFTracking/interface/PFTrackAlgoTools.h"
0049 #include "RecoParticleFlow/PFTracking/interface/PFTrackTransformer.h"
0050 #include "TrackingTools/GsfTools/interface/MultiTrajectoryStateMode.h"
0051 #include "TrackingTools/GsfTools/interface/MultiTrajectoryStateTransform.h"
0052 #include "TrackingTools/PatternTools/interface/Trajectory.h"
0053 #include "TrackingTools/Records/interface/TransientTrackRecord.h"
0054 #include "TrackingTools/TransientTrack/interface/TransientTrackBuilder.h"
0055
0056 namespace {
0057
0058   constexpr float square(float x) { return x * x; };
0059
0060 }  // namespace
0061
0062 class PFElecTkProducer final : public edm::stream::EDProducer<edm::GlobalCache<convbremhelpers::HeavyObjectCache> > {
0063 public:
0064   ///Constructor
0065   explicit PFElecTkProducer(const edm::ParameterSet&, const convbremhelpers::HeavyObjectCache*);
0066
0067   static std::unique_ptr<convbremhelpers::HeavyObjectCache> initializeGlobalCache(const edm::ParameterSet& conf) {
0068     return std::make_unique<convbremhelpers::HeavyObjectCache>(conf);
0069   }
0070
0071   static void globalEndJob(convbremhelpers::HeavyObjectCache const*) {}
0072
0073 private:
0074   void beginRun(const edm::Run&, const edm::EventSetup&) override;
0075   void endRun(const edm::Run&, const edm::EventSetup&) override;
0076
0077   ///Produce the PFRecTrack collection
0078   void produce(edm::Event&, const edm::EventSetup&) override;
0079
0080   int findPfRef(const reco::PFRecTrackCollection& pfRTkColl,
0081                 const reco::GsfTrack&,
0082                 edm::soa::EtaPhiTableView trackEtaPhiTable);
0083
0084   bool applySelection(const reco::GsfTrack&);
0085
0086   bool resolveGsfTracks(const std::vector<reco::GsfPFRecTrack>& GsfPFVec,
0087                         unsigned int ngsf,
0088                         std::vector<unsigned int>& secondaries,
0089                         const reco::PFClusterCollection& theEClus);
0090
0091   float minTangDist(const reco::GsfPFRecTrack& primGsf, const reco::GsfPFRecTrack& secGsf);
0092
0093   bool isSameEgSC(const reco::ElectronSeed& nSeed,
0094                   const reco::ElectronSeed& iSeed,
0095                   bool& bothGsfEcalDriven,
0096                   float& SCEnergy);
0097
0098   bool isSharingEcalEnergyWithEgSC(const reco::GsfPFRecTrack& nGsfPFRecTrack,
0099                                    const reco::GsfPFRecTrack& iGsfPFRecTrack,
0100                                    const reco::ElectronSeed& nSeed,
0101                                    const reco::ElectronSeed& iSeed,
0102                                    const reco::PFClusterCollection& theEClus,
0103                                    bool& bothGsfTrackerDriven,
0104                                    bool& nEcalDriven,
0105                                    bool& iEcalDriven,
0106                                    float& nEnergy,
0107                                    float& iEnergy);
0108
0109   bool isInnerMost(const reco::GsfTrackRef& nGsfTrack, const reco::GsfTrackRef& iGsfTrack, bool& sameLayer);
0110
0111   bool isInnerMostWithLostHits(const reco::GsfTrackRef& nGsfTrack, const reco::GsfTrackRef& iGsfTrack, bool& sameLayer);
0112
0113   void createGsfPFRecTrackRef(const edm::OrphanHandle<reco::GsfPFRecTrackCollection>& gsfPfHandle,
0114                               std::vector<reco::GsfPFRecTrack>& gsfPFRecTrackPrimary,
0115                               const std::map<unsigned int, std::vector<reco::GsfPFRecTrack> >& MapPrimSec);
0116
0117   // ----------member data ---------------------------
0118   reco::GsfPFRecTrack pftrack_;
0119   reco::GsfPFRecTrack secpftrack_;
0120   edm::ParameterSet conf_;
0121   edm::EDGetTokenT<reco::GsfTrackCollection> gsfTrackLabel_;
0122   edm::EDGetTokenT<reco::PFRecTrackCollection> pfTrackLabel_;
0123   edm::EDGetTokenT<reco::VertexCollection> primVtxLabel_;
0124   edm::EDGetTokenT<reco::PFClusterCollection> pfEcalClusters_;
0125   edm::EDGetTokenT<reco::PFDisplacedTrackerVertexCollection> pfNuclear_;
0126   edm::EDGetTokenT<reco::PFConversionCollection> pfConv_;
0127   edm::EDGetTokenT<reco::PFV0Collection> pfV0_;
0128   bool useNuclear_;
0129   bool useConversions_;
0130   bool useV0_;
0131   bool applyAngularGsfClean_;
0132   double detaCutGsfClean_;
0133   double dphiCutGsfClean_;
0134
0135   const edm::ESGetToken<MagneticField, IdealMagneticFieldRecord> magFieldToken_;
0136   const edm::ESGetToken<TrackerGeometry, TrackerDigiGeometryRecord> tkerGeomToken_;
0137   const edm::ESGetToken<TransientTrackBuilder, TransientTrackRecord> transientTrackToken_;
0138
0139   ///PFTrackTransformer
0140   std::unique_ptr<PFTrackTransformer> pfTransformer_;
0141   MultiTrajectoryStateTransform mtsTransform_;
0142   std::unique_ptr<ConvBremPFTrackFinder> convBremFinder_;
0143
0144   ///Trajectory of GSfTracks in the event?
0145   bool trajinev_;
0146   bool modemomentum_;
0147   bool applySel_;
0148   bool applyGsfClean_;
0149   bool useFifthStepForEcalDriven_;
0150   bool useFifthStepForTrackDriven_;
0151   //   bool useFifthStepSec_;
0152   bool debugGsfCleaning_;
0153   double SCEne_;
0154   double detaGsfSC_;
0155   double dphiGsfSC_;
0156   double maxPtConvReco_;
0157
0158   /// Conv Brem Finder
0159   bool useConvBremFinder_;
0160
0161   double mvaConvBremFinderIDBarrelLowPt_;
0162   double mvaConvBremFinderIDBarrelHighPt_;
0163   double mvaConvBremFinderIDEndcapsLowPt_;
0164   double mvaConvBremFinderIDEndcapsHighPt_;
0165   std::string path_mvaWeightFileConvBremBarrelLowPt_;
0166   std::string path_mvaWeightFileConvBremBarrelHighPt_;
0167   std::string path_mvaWeightFileConvBremEndcapsLowPt_;
0168   std::string path_mvaWeightFileConvBremEndcapsHighPt_;
0169
0170   // cache for multitrajectory states
0171   std::vector<double> gsfInnerMomentumCache_;
0172 };
0173
0174 using namespace std;
0175 using namespace edm;
0176 using namespace reco;
0177
0178 PFElecTkProducer::PFElecTkProducer(const ParameterSet& iConfig, const convbremhelpers::HeavyObjectCache*)
0179     : conf_(iConfig),
0180       magFieldToken_(esConsumes<edm::Transition::BeginRun>()),
0181       tkerGeomToken_(esConsumes<edm::Transition::BeginRun>()),
0182       transientTrackToken_(esConsumes<edm::Transition::BeginRun>(edm::ESInputTag("", "TransientTrackBuilder"))) {
0183   gsfTrackLabel_ = consumes<reco::GsfTrackCollection>(iConfig.getParameter<InputTag>("GsfTrackModuleLabel"));
0184
0185   pfTrackLabel_ = consumes<reco::PFRecTrackCollection>(iConfig.getParameter<InputTag>("PFRecTrackLabel"));
0186
0187   primVtxLabel_ = consumes<reco::VertexCollection>(iConfig.getParameter<InputTag>("PrimaryVertexLabel"));
0188
0189   pfEcalClusters_ = consumes<reco::PFClusterCollection>(iConfig.getParameter<InputTag>("PFEcalClusters"));
0190
0191   pfNuclear_ = consumes<reco::PFDisplacedTrackerVertexCollection>(iConfig.getParameter<InputTag>("PFNuclear"));
0192
0193   pfConv_ = consumes<reco::PFConversionCollection>(iConfig.getParameter<InputTag>("PFConversions"));
0194
0195   pfV0_ = consumes<reco::PFV0Collection>(iConfig.getParameter<InputTag>("PFV0"));
0196
0197   useNuclear_ = iConfig.getParameter<bool>("useNuclear");
0198   useConversions_ = iConfig.getParameter<bool>("useConversions");
0199   useV0_ = iConfig.getParameter<bool>("useV0");
0200   debugGsfCleaning_ = iConfig.getParameter<bool>("debugGsfCleaning");
0201
0202   produces<GsfPFRecTrackCollection>();
0203   produces<GsfPFRecTrackCollection>("Secondary").setBranchAlias("secondary");
0204
0205   trajinev_ = iConfig.getParameter<bool>("TrajInEvents");
0206   modemomentum_ = iConfig.getParameter<bool>("ModeMomentum");
0207   applySel_ = iConfig.getParameter<bool>("applyEGSelection");
0208   applyGsfClean_ = iConfig.getParameter<bool>("applyGsfTrackCleaning");
0209   applyAngularGsfClean_ = iConfig.getParameter<bool>("applyAlsoGsfAngularCleaning");
0210   detaCutGsfClean_ = iConfig.getParameter<double>("maxDEtaGsfAngularCleaning");
0211   dphiCutGsfClean_ = iConfig.getParameter<double>("maxDPhiBremTangGsfAngularCleaning");
0212   useFifthStepForTrackDriven_ = iConfig.getParameter<bool>("useFifthStepForTrackerDrivenGsf");
0213   useFifthStepForEcalDriven_ = iConfig.getParameter<bool>("useFifthStepForEcalDrivenGsf");
0214   maxPtConvReco_ = iConfig.getParameter<double>("MaxConvBremRecoPT");
0215   detaGsfSC_ = iConfig.getParameter<double>("MinDEtaGsfSC");
0216   dphiGsfSC_ = iConfig.getParameter<double>("MinDPhiGsfSC");
0217   SCEne_ = iConfig.getParameter<double>("MinSCEnergy");
0218
0219   // set parameter for convBremFinder
0220   useConvBremFinder_ = iConfig.getParameter<bool>("useConvBremFinder");
0221
0222   mvaConvBremFinderIDBarrelLowPt_ = iConfig.getParameter<double>("pf_convBremFinderID_mvaCutBarrelLowPt");
0223   mvaConvBremFinderIDBarrelHighPt_ = iConfig.getParameter<double>("pf_convBremFinderID_mvaCutBarrelHighPt");
0224   mvaConvBremFinderIDEndcapsLowPt_ = iConfig.getParameter<double>("pf_convBremFinderID_mvaCutEndcapsLowPt");
0225   mvaConvBremFinderIDEndcapsHighPt_ = iConfig.getParameter<double>("pf_convBremFinderID_mvaCutEndcapsHighPt");
0226 }
0227
0228 //
0229 // member functions
0230 //
0231
0232 // ------------ method called to produce the data  ------------
0233 void PFElecTkProducer::produce(Event& iEvent, const EventSetup& iSetup) {
0234   //create the empty collections
0235   auto gsfPFRecTrackCollection = std::make_unique<GsfPFRecTrackCollection>();
0236
0237   auto gsfPFRecTrackCollectionSecondary = std::make_unique<GsfPFRecTrackCollection>();
0238
0239   //read collections of tracks
0240   Handle<GsfTrackCollection> gsftrackscoll;
0241   iEvent.getByToken(gsfTrackLabel_, gsftrackscoll);
0242
0243   //read collections of trajectories
0244   Handle<vector<Trajectory> > TrajectoryCollection;
0245
0246   //read pfrectrack collection
0247   Handle<PFRecTrackCollection> thePfRecTrackCollection;
0248   iEvent.getByToken(pfTrackLabel_, thePfRecTrackCollection);
0249
0250   // SoA structure for frequently used track information.
0251   // LazyConstructed so it is only filled when needed, i.e., when there is an electron in the event.
0252   auto trackEtaPhiTable = makeLazy<edm::soa::EtaPhiTable>(*thePfRecTrackCollection);
0253
0254   // PFClusters
0255   Handle<PFClusterCollection> theECPfClustCollection;
0256   iEvent.getByToken(pfEcalClusters_, theECPfClustCollection);
0257   const PFClusterCollection& theEcalClusters = *(theECPfClustCollection.product());
0258
0259   //Primary Vertexes
0260   Handle<reco::VertexCollection> thePrimaryVertexColl;
0261   iEvent.getByToken(primVtxLabel_, thePrimaryVertexColl);
0262
0263   // Displaced Vertex
0264   Handle<reco::PFDisplacedTrackerVertexCollection> pfNuclears;
0265   if (useNuclear_)
0266     iEvent.getByToken(pfNuclear_, pfNuclears);
0267
0268   // Conversions
0269   Handle<reco::PFConversionCollection> pfConversions;
0270   if (useConversions_)
0271     iEvent.getByToken(pfConv_, pfConversions);
0272
0273   // V0
0274   Handle<reco::PFV0Collection> pfV0;
0275   if (useV0_)
0276     iEvent.getByToken(pfV0_, pfV0);
0277
0278   GsfTrackCollection gsftracks = *(gsftrackscoll.product());
0279   vector<Trajectory> tjvec(0);
0280   if (trajinev_) {
0281     iEvent.getByToken(gsfTrackLabel_, TrajectoryCollection);
0282
0283     tjvec = *(TrajectoryCollection.product());
0284   }
0285
0286   vector<reco::GsfPFRecTrack> selGsfPFRecTracks;
0287   vector<reco::GsfPFRecTrack> primaryGsfPFRecTracks;
0288   std::map<unsigned int, std::vector<reco::GsfPFRecTrack> > GsfPFMap;
0289
0290   for (unsigned igsf = 0; igsf < gsftracks.size(); igsf++) {
0291     GsfTrackRef trackRef(gsftrackscoll, igsf);
0292     gsfInnerMomentumCache_.push_back(trackRef->pMode());
0293     TrajectoryStateOnSurface i_inTSOS = mtsTransform_.innerStateOnSurface((*trackRef));
0294     GlobalVector i_innMom;
0295     if (i_inTSOS.isValid()) {
0296       multiTrajectoryStateMode::momentumFromModeCartesian(i_inTSOS, i_innMom);
0297       gsfInnerMomentumCache_.back() = i_innMom.mag();
0298     }
0299   }
0300
0301   for (unsigned int igsf = 0; igsf < gsftracks.size(); igsf++) {
0302     GsfTrackRef trackRef(gsftrackscoll, igsf);
0303
0304     int kf_ind = findPfRef(*thePfRecTrackCollection, gsftracks[igsf], trackEtaPhiTable.value());
0305
0306     if (kf_ind >= 0) {
0307       PFRecTrackRef kf_ref(thePfRecTrackCollection, kf_ind);
0308
0309       // remove fifth step tracks
0310       if (useFifthStepForEcalDriven_ == false || useFifthStepForTrackDriven_ == false) {
0311         bool isFifthStepTrack = PFTrackAlgoTools::isFifthStep(kf_ref->trackRef()->algo());
0312         bool isEcalDriven = true;
0313         bool isTrackerDriven = true;
0314
0315         if (trackRef->seedRef().get() == nullptr) {
0316           isEcalDriven = false;
0317           isTrackerDriven = false;
0318         } else {
0319           auto const& SeedFromRef = static_cast<ElectronSeed const&>(*(trackRef->extra()->seedRef()));
0320           if (SeedFromRef.caloCluster().isNull())
0321             isEcalDriven = false;
0322           if (SeedFromRef.ctfTrack().isNull())
0323             isTrackerDriven = false;
0324         }
0325         //note: the same track could be both ecalDriven and trackerDriven
0326         if (isFifthStepTrack && isEcalDriven && isTrackerDriven == false && useFifthStepForEcalDriven_ == false) {
0327           continue;
0328         }
0329
0330         if (isFifthStepTrack && isTrackerDriven && isEcalDriven == false && useFifthStepForTrackDriven_ == false) {
0331           continue;
0332         }
0333
0334         if (isFifthStepTrack && isTrackerDriven && isEcalDriven && useFifthStepForTrackDriven_ == false &&
0335             useFifthStepForEcalDriven_ == false) {
0336           continue;
0337         }
0338       }
0339
0340       pftrack_ = GsfPFRecTrack(gsftracks[igsf].charge(), reco::PFRecTrack::GSF, igsf, trackRef, kf_ref);
0341     } else {
0342       PFRecTrackRef dummyRef;
0343       pftrack_ = GsfPFRecTrack(gsftracks[igsf].charge(), reco::PFRecTrack::GSF, igsf, trackRef, dummyRef);
0344     }
0345
0346     bool validgsfbrem = false;
0347     if (trajinev_) {
0348       validgsfbrem = pfTransformer_->addPointsAndBrems(pftrack_, gsftracks[igsf], tjvec[igsf], modemomentum_);
0349     } else {
0350       validgsfbrem = pfTransformer_->addPointsAndBrems(pftrack_, gsftracks[igsf], mtsTransform_);
0351     }
0352
0353     bool passSel = true;
0354     if (applySel_)
0355       passSel = applySelection(gsftracks[igsf]);
0356
0357     if (validgsfbrem && passSel)
0358       selGsfPFRecTracks.push_back(pftrack_);
0359   }
0360
0361   unsigned int count_primary = 0;
0362   if (!selGsfPFRecTracks.empty()) {
0363     for (unsigned int ipfgsf = 0; ipfgsf < selGsfPFRecTracks.size(); ipfgsf++) {
0364       vector<unsigned int> secondaries(0);
0365       secondaries.clear();
0366       bool keepGsf = true;
0367
0368       if (applyGsfClean_) {
0369         keepGsf = resolveGsfTracks(selGsfPFRecTracks, ipfgsf, secondaries, theEcalClusters);
0370       }
0371
0372       //is primary?
0373       if (keepGsf == true) {
0374         // Find kf tracks from converted brem photons
0375         if (convBremFinder_->foundConvBremPFRecTrack(thePfRecTrackCollection,
0376                                                      thePrimaryVertexColl,
0377                                                      pfNuclears,
0378                                                      pfConversions,
0379                                                      pfV0,
0380                                                      globalCache(),
0381                                                      useNuclear_,
0382                                                      useConversions_,
0383                                                      useV0_,
0384                                                      theEcalClusters,
0385                                                      selGsfPFRecTracks[ipfgsf])) {
0386           const vector<PFRecTrackRef>& convBremPFRecTracks(convBremFinder_->getConvBremPFRecTracks());
0387           for (unsigned int ii = 0; ii < convBremPFRecTracks.size(); ii++) {
0388             selGsfPFRecTracks[ipfgsf].addConvBremPFRecTrackRef(convBremPFRecTracks[ii]);
0389           }
0390         }
0391
0392         // save primaries gsf tracks
0393         //  gsfPFRecTrackCollection->push_back(selGsfPFRecTracks[ipfgsf]);
0394         primaryGsfPFRecTracks.push_back(selGsfPFRecTracks[ipfgsf]);
0395
0396         // NOTE:: THE TRACKID IS USED TO LINK THE PRIMARY GSF TRACK. THIS NEEDS
0397         // TO BE CHANGED AS SOON AS IT IS POSSIBLE TO CHANGE DATAFORMATS
0398         // A MODIFICATION HERE IMPLIES A MODIFICATION IN PFBLOCKALGO.CC/H
0399         unsigned int primGsfIndex = selGsfPFRecTracks[ipfgsf].trackId();
0400         vector<reco::GsfPFRecTrack> trueGsfPFRecTracks;
0401         if (!secondaries.empty()) {
0402           // loop on secondaries gsf tracks (from converted brems)
0403           for (unsigned int isecpfgsf = 0; isecpfgsf < secondaries.size(); isecpfgsf++) {
0404             PFRecTrackRef refsecKF = selGsfPFRecTracks[(secondaries[isecpfgsf])].kfPFRecTrackRef();
0405
0406             unsigned int secGsfIndex = selGsfPFRecTracks[(secondaries[isecpfgsf])].trackId();
0407             GsfTrackRef secGsfRef = selGsfPFRecTracks[(secondaries[isecpfgsf])].gsfTrackRef();
0408
0409             if (refsecKF.isNonnull()) {
0410               // NOTE::IT SAVED THE TRACKID OF THE PRIMARY!!! THIS IS USED IN PFBLOCKALGO.CC/H
0411               secpftrack_ = GsfPFRecTrack(
0412                   gsftracks[secGsfIndex].charge(), reco::PFRecTrack::GSF, primGsfIndex, secGsfRef, refsecKF);
0413             } else {
0414               PFRecTrackRef dummyRef;
0415               // NOTE::IT SAVED THE TRACKID OF THE PRIMARY!!! THIS IS USED IN PFBLOCKALGO.CC/H
0416               secpftrack_ = GsfPFRecTrack(
0417                   gsftracks[secGsfIndex].charge(), reco::PFRecTrack::GSF, primGsfIndex, secGsfRef, dummyRef);
0418             }
0419
0420             bool validgsfbrem = false;
0421             if (trajinev_) {
0422               validgsfbrem = pfTransformer_->addPointsAndBrems(
0423                   secpftrack_, gsftracks[secGsfIndex], tjvec[secGsfIndex], modemomentum_);
0424             } else {
0425               validgsfbrem = pfTransformer_->addPointsAndBrems(secpftrack_, gsftracks[secGsfIndex], mtsTransform_);
0426             }
0427
0428             if (validgsfbrem) {
0429               gsfPFRecTrackCollectionSecondary->push_back(secpftrack_);
0430               trueGsfPFRecTracks.push_back(secpftrack_);
0431             }
0432           }
0433         }
0434         GsfPFMap.insert(pair<unsigned int, std::vector<reco::GsfPFRecTrack> >(count_primary, trueGsfPFRecTracks));
0435         trueGsfPFRecTracks.clear();
0436         count_primary++;
0437       }
0438     }
0439   }
0440
0441   const edm::OrphanHandle<GsfPFRecTrackCollection> gsfPfRefProd =
0442       iEvent.put(std::move(gsfPFRecTrackCollectionSecondary), "Secondary");
0443
0444   //now the secondary GsfPFRecTracks are in the event, the Ref can be created
0445   createGsfPFRecTrackRef(gsfPfRefProd, primaryGsfPFRecTracks, GsfPFMap);
0446
0447   for (unsigned int iGSF = 0; iGSF < primaryGsfPFRecTracks.size(); iGSF++) {
0448     gsfPFRecTrackCollection->push_back(primaryGsfPFRecTracks[iGSF]);
0449   }
0450   iEvent.put(std::move(gsfPFRecTrackCollection));
0451
0452   selGsfPFRecTracks.clear();
0453   GsfPFMap.clear();
0454   primaryGsfPFRecTracks.clear();
0455
0456   std::vector<double>().swap(gsfInnerMomentumCache_);
0457 }
0458
0459 // create the secondary GsfPFRecTracks Ref
0460 void PFElecTkProducer::createGsfPFRecTrackRef(
0461     const edm::OrphanHandle<reco::GsfPFRecTrackCollection>& gsfPfHandle,
0462     std::vector<reco::GsfPFRecTrack>& gsfPFRecTrackPrimary,
0463     const std::map<unsigned int, std::vector<reco::GsfPFRecTrack> >& MapPrimSec) {
0464   unsigned int cgsf = 0;
0465   unsigned int csecgsf = 0;
0466   for (std::map<unsigned int, std::vector<reco::GsfPFRecTrack> >::const_iterator igsf = MapPrimSec.begin();
0467        igsf != MapPrimSec.end();
0468        igsf++, cgsf++) {
0469     vector<reco::GsfPFRecTrack> SecGsfPF = igsf->second;
0470     for (unsigned int iSecGsf = 0; iSecGsf < SecGsfPF.size(); iSecGsf++) {
0471       edm::Ref<reco::GsfPFRecTrackCollection> refgprt(gsfPfHandle, csecgsf);
0472       gsfPFRecTrackPrimary[cgsf].addConvBremGsfPFRecTrackRef(refgprt);
0473       ++csecgsf;
0474     }
0475   }
0476
0477   return;
0478 }
0479 // ------------- method for find the corresponding kf pfrectrack ---------------------
0480 int PFElecTkProducer::findPfRef(const reco::PFRecTrackCollection& pfRTkColl,
0481                                 const reco::GsfTrack& gsftk,
0482                                 edm::soa::EtaPhiTableView trackEtaPhiTable) {
0483   if (gsftk.seedRef().get() == nullptr) {
0484     return -1;
0485   }
0486
0487   // maximum delta_r2 for matching
0488   constexpr float maxDR2 = square(0.05f);
0489
0490   // precompute expensive trigonometry
0491   float gsftkEta = gsftk.eta();
0492   float gsftkPhi = gsftk.phi();
0493   auto const& gsftkHits = gsftk.seedRef()->recHits();
0494
0495   auto const& electronSeedFromRef = static_cast<ElectronSeed const&>(*(gsftk.extra()->seedRef()));
0496   //CASE 1 ELECTRONSEED DOES NOT HAVE A REF TO THE CKFTRACK
0497   if (electronSeedFromRef.ctfTrack().isNull()) {
0498     unsigned int i_pf = 0;
0499     int ibest = -1;
0500     unsigned int ish_max = 0;
0501     float dr2_min = square(1000.f);
0502
0503     //SEARCH THE PFRECTRACK THAT SHARES HITS WITH THE ELECTRON SEED
0504     // Here the cpu time has been be improved.
0505     for (auto const& pft : trackEtaPhiTable) {
0506       unsigned int ish = 0;
0507
0508       using namespace edm::soa::col;
0509       const float dr2 = reco::deltaR2(pft.get<Eta>(), pft.get<Phi>(), gsftkEta, gsftkPhi);
0510
0511       if (dr2 <= maxDR2) {
0512         for (auto const& hhit : pfRTkColl[i_pf].trackRef()->recHits()) {
0513           if (!hhit->isValid())
0514             continue;
0515           for (auto const& hit : gsftkHits) {
0516             if (hit.isValid() && hhit->sharesInput(&hit, TrackingRecHit::all))
0517               ish++;
0518           }
0519         }
0520
0521         if ((ish > ish_max) || ((ish == ish_max) && (dr2 < dr2_min))) {
0522           ish_max = ish;
0523           dr2_min = dr2;
0524           ibest = i_pf;
0525         }
0526       }
0527
0528       i_pf++;
0529     }
0530
0531     return ((ish_max == 0) || (dr2_min > maxDR2)) ? -1 : ibest;
0532   } else {
0533     //ELECTRON SEED HAS A REFERENCE
0534
0535     unsigned int i_pf = 0;
0536
0537     for (auto const& pft : pfRTkColl) {
0538       //REF COMPARISON
0539       if (pft.trackRef() == electronSeedFromRef.ctfTrack()) {
0540         return i_pf;
0541       }
0542       i_pf++;
0543     }
0544   }
0545   return -1;
0546 }
0547
0548 // -- method to apply gsf electron selection to EcalDriven seeds
0549 bool PFElecTkProducer::applySelection(const reco::GsfTrack& gsftk) {
0550   if (gsftk.seedRef().get() == nullptr)
0551     return false;
0552   auto const& ElSeedFromRef = static_cast<ElectronSeed const&>(*(gsftk.extra()->seedRef()));
0553
0554   bool passCut = false;
0555   if (ElSeedFromRef.ctfTrack().isNull()) {
0556     if (ElSeedFromRef.caloCluster().isNull())
0557       return passCut;
0558     auto const* scRef = static_cast<SuperCluster const*>(ElSeedFromRef.caloCluster().get());
0559     //do this just to know if exist a SC?
0560     if (scRef) {
0561       float caloEne = scRef->energy();
0562       float feta = fabs(scRef->eta() - gsftk.etaMode());
0563       float fphi = fabs(scRef->phi() - gsftk.phiMode());
0564       if (fphi > TMath::Pi())
0565         fphi -= TMath::TwoPi();
0566       if (caloEne > SCEne_ && feta < detaGsfSC_ && fabs(fphi) < dphiGsfSC_)
0567         passCut = true;
0568     }
0569   } else {
0570     // get all the gsf found by tracker driven
0571     passCut = true;
0572   }
0573   return passCut;
0574 }
0575 bool PFElecTkProducer::resolveGsfTracks(const vector<reco::GsfPFRecTrack>& GsfPFVec,
0576                                         unsigned int ngsf,
0577                                         vector<unsigned int>& secondaries,
0578                                         const reco::PFClusterCollection& theEClus) {
0579   bool debugCleaning = debugGsfCleaning_;
0580   bool n_keepGsf = true;
0581
0582   const reco::GsfTrackRef& nGsfTrack = GsfPFVec[ngsf].gsfTrackRef();
0583
0584   if (nGsfTrack->seedRef().get() == nullptr)
0585     return false;
0586   auto const& nElSeedFromRef = static_cast<ElectronSeed const&>(*(nGsfTrack->extra()->seedRef()));
0587
0588   /* // now gotten from cache below
0589   TrajectoryStateOnSurface inTSOS = mtsTransform_.innerStateOnSurface((*nGsfTrack));
0590   GlobalVector ninnMom;
0591   float nPin =  nGsfTrack->pMode();
0592   if(inTSOS.isValid()){
0593     multiTrajectoryStateMode::momentumFromModeCartesian(inTSOS,ninnMom);
0594     nPin = ninnMom.mag();
0595   }
0596   */
0597   float nPin = gsfInnerMomentumCache_[nGsfTrack.key()];
0598
0599   float neta = nGsfTrack->innerMomentum().eta();
0600   float nphi = nGsfTrack->innerMomentum().phi();
0601
0602   if (debugCleaning)
0603     cout << " PFElecTkProducer:: considering track " << nGsfTrack->pt() << " eta,phi " << nGsfTrack->eta() << ", "
0604          << nGsfTrack->phi() << endl;
0605
0606   for (unsigned int igsf = 0; igsf < GsfPFVec.size(); igsf++) {
0607     if (igsf != ngsf) {
0608       reco::GsfTrackRef iGsfTrack = GsfPFVec[igsf].gsfTrackRef();
0609
0610       if (debugCleaning)
0611         cout << " PFElecTkProducer:: and  comparing with track " << iGsfTrack->pt() << " eta,phi " << iGsfTrack->eta()
0612              << ", " << iGsfTrack->phi() << endl;
0613
0614       float ieta = iGsfTrack->innerMomentum().eta();
0615       float iphi = iGsfTrack->innerMomentum().phi();
0616       float feta = fabs(neta - ieta);
0617       float fphi = fabs(nphi - iphi);
0618       if (fphi > TMath::Pi())
0619         fphi -= TMath::TwoPi();
0620
0621       // apply a superloose preselection only to avoid un-useful cpu time: hard-coded for this reason
0622       if (feta < 0.5 && fabs(fphi) < 1.0) {
0623         if (debugCleaning)
0624           cout << " Entering angular superloose preselection " << endl;
0625
0626         /* //now taken from cache below
0627         TrajectoryStateOnSurface i_inTSOS = mtsTransform_.innerStateOnSurface((*iGsfTrack));
0628         GlobalVector i_innMom;
0629         float iPin = iGsfTrack->pMode();
0630         if(i_inTSOS.isValid()){
0631           multiTrajectoryStateMode::momentumFromModeCartesian(i_inTSOS,i_innMom);
0632           iPin = i_innMom.mag();
0633         }
0634         */
0635         float iPin = gsfInnerMomentumCache_[iGsfTrack.key()];
0636
0637         if (iGsfTrack->seedRef().get() == nullptr)
0638           continue;
0639         auto const& iElSeedFromRef = static_cast<ElectronSeed const&>(*(iGsfTrack->extra()->seedRef()));
0640
0641         float SCEnergy = -1.;
0642         // Check if two tracks match the same SC
0643         bool areBothGsfEcalDriven = false;
0644         ;
0645         bool isSameSC = isSameEgSC(nElSeedFromRef, iElSeedFromRef, areBothGsfEcalDriven, SCEnergy);
0646
0647         // CASE1 both GsfTracks ecalDriven and match the same SC
0648         if (areBothGsfEcalDriven) {
0649           if (isSameSC) {
0650             float nEP = SCEnergy / nPin;
0651             float iEP = SCEnergy / iPin;
0652             if (debugCleaning)
0653               cout << " Entering SAME supercluster case "
0654                    << " nEP " << nEP << " iEP " << iEP << endl;
0655
0656             // if same SC take the closest or if same
0657             // distance the best E/p
0658
0659             // Innermost using LostHits technology
0660             bool isSameLayer = false;
0661             bool iGsfInnermostWithLostHits = isInnerMostWithLostHits(nGsfTrack, iGsfTrack, isSameLayer);
0662
0663             if (debugCleaning)
0664               cout << " iGsf is InnerMostWithLostHits " << iGsfInnermostWithLostHits << " isSameLayer " << isSameLayer
0665                    << endl;
0666
0667             if (iGsfInnermostWithLostHits) {
0668               n_keepGsf = false;
0669               return n_keepGsf;
0670             } else if (isSameLayer) {
0671               if (fabs(iEP - 1) < fabs(nEP - 1)) {
0672                 n_keepGsf = false;
0673                 return n_keepGsf;
0674               } else {
0675                 secondaries.push_back(igsf);
0676               }
0677             } else {
0678               // save secondaries gsf track (put selection)
0679               secondaries.push_back(igsf);
0680             }
0681           }  // end same SC case
0682         } else {
0683           // enter in the condition where at least one track is trackerDriven
0684           float minBremDphi = minTangDist(GsfPFVec[ngsf], GsfPFVec[igsf]);
0685           float nETot = 0.;
0686           float iETot = 0.;
0687           bool isBothGsfTrackerDriven = false;
0688           bool nEcalDriven = false;
0689           bool iEcalDriven = false;
0690           bool isSameScEgPf = isSharingEcalEnergyWithEgSC(GsfPFVec[ngsf],
0691                                                           GsfPFVec[igsf],
0692                                                           nElSeedFromRef,
0693                                                           iElSeedFromRef,
0694                                                           theEClus,
0695                                                           isBothGsfTrackerDriven,
0696                                                           nEcalDriven,
0697                                                           iEcalDriven,
0698                                                           nETot,
0699                                                           iETot);
0700
0701           // check if the first hit of iGsfTrack < nGsfTrack
0702           bool isSameLayer = false;
0703           bool iGsfInnermostWithLostHits = isInnerMostWithLostHits(nGsfTrack, iGsfTrack, isSameLayer);
0704
0705           if (isSameScEgPf) {
0706             // CASE 2 : One Gsf has reference to a SC and the other one not or both not
0707
0708             if (debugCleaning) {
0709               cout << " Sharing ECAL energy passed "
0710                    << " nEtot " << nETot << " iEtot " << iETot << endl;
0711               if (isBothGsfTrackerDriven)
0712                 cout << " Both Track are trackerDriven " << endl;
0713             }
0714
0715             // Innermost using LostHits technology
0716             if (iGsfInnermostWithLostHits) {
0717               n_keepGsf = false;
0718               return n_keepGsf;
0719             } else if (isSameLayer) {
0720               // Thirt Case:  One Gsf has reference to a SC and the other one not or both not
0721               // gsf tracks starts from the same layer
0722               // check number of sharing modules (at least 50%)
0723               // check number of sharing hits (at least 2)
0724               // check charge flip inner/outer
0725
0726               // they share energy
0727               if (isBothGsfTrackerDriven == false) {
0728                 // if at least one Gsf track is EcalDriven choose that one.
0729                 if (iEcalDriven) {
0730                   n_keepGsf = false;
0731                   return n_keepGsf;
0732                 } else {
0733                   secondaries.push_back(igsf);
0734                 }
0735               } else {
0736                 // if both tracks are tracker driven choose that one with the best E/p
0737                 // with ETot = max(En,Ei)
0738
0739                 float ETot = -1;
0740                 if (nETot != iETot) {
0741                   if (nETot > iETot)
0742                     ETot = nETot;
0743                   else
0744                     ETot = iETot;
0745                 } else {
0746                   ETot = nETot;
0747                 }
0748                 float nEP = ETot / nPin;
0749                 float iEP = ETot / iPin;
0750
0751                 if (debugCleaning)
0752                   cout << " nETot " << nETot << " iETot " << iETot << " ETot " << ETot << endl
0753                        << " nPin " << nPin << " iPin " << iPin << " nEP " << nEP << " iEP " << iEP << endl;
0754
0755                 if (fabs(iEP - 1) < fabs(nEP - 1)) {
0756                   n_keepGsf = false;
0757                   return n_keepGsf;
0758                 } else {
0759                   secondaries.push_back(igsf);
0760                 }
0761               }
0762             } else {
0763               secondaries.push_back(igsf);
0764             }
0765           } else if (feta < detaCutGsfClean_ && minBremDphi < dphiCutGsfClean_) {
0766             // very close tracks
0767             bool secPushedBack = false;
0768             if (nEcalDriven == false && nETot == 0.) {
0769               n_keepGsf = false;
0770               return n_keepGsf;
0771             } else if (iEcalDriven == false && iETot == 0.) {
0772               secondaries.push_back(igsf);
0773               secPushedBack = true;
0774             }
0775             if (debugCleaning)
0776               cout << " Close Tracks "
0777                    << " feta " << feta << " fabs(fphi) " << fabs(fphi) << " minBremDphi " << minBremDphi << " nETot "
0778                    << nETot << " iETot " << iETot << " nLostHits " << nGsfTrack->missingInnerHits() << " iLostHits "
0779                    << iGsfTrack->missingInnerHits() << endl;
0780
0781             // apply selection only if one track has lost hits
0782             if (applyAngularGsfClean_) {
0783               if (iGsfInnermostWithLostHits) {
0784                 n_keepGsf = false;
0785                 return n_keepGsf;
0786               } else if (isSameLayer == false) {
0787                 if (secPushedBack == false)
0788                   secondaries.push_back(igsf);
0789               }
0790             }
0791           } else if (feta < 0.1 && minBremDphi < 0.2) {
0792             // failed all the conditions, discard only tracker driven tracks
0793             // with no PFClusters linked.
0794             if (debugCleaning)
0795               cout << " Close Tracks and failed all the conditions "
0796                    << " feta " << feta << " fabs(fphi) " << fabs(fphi) << " minBremDphi " << minBremDphi << " nETot "
0797                    << nETot << " iETot " << iETot << " nLostHits " << nGsfTrack->missingInnerHits() << " iLostHits "
0798                    << iGsfTrack->missingInnerHits() << endl;
0799
0800             if (nEcalDriven == false && nETot == 0.) {
0801               n_keepGsf = false;
0802               return n_keepGsf;
0803             }
0804             // Here I do not push back the secondary because considered fakes...
0805           }
0806         }
0807       }
0808     }
0809   }
0810
0811   return n_keepGsf;
0812 }
0813 float PFElecTkProducer::minTangDist(const reco::GsfPFRecTrack& primGsf, const reco::GsfPFRecTrack& secGsf) {
0814   float minDphi = 1000.;
0815
0816   std::vector<reco::PFBrem> primPFBrem = primGsf.PFRecBrem();
0817   std::vector<reco::PFBrem> secPFBrem = secGsf.PFRecBrem();
0818
0819   unsigned int cbrem = 0;
0820   for (unsigned isbrem = 0; isbrem < secPFBrem.size(); isbrem++) {
0821     if (secPFBrem[isbrem].indTrajPoint() == 99)
0822       continue;
0823     const reco::PFTrajectoryPoint& atSecECAL =
0824         secPFBrem[isbrem].extrapolatedPoint(reco::PFTrajectoryPoint::ECALEntrance);
0825     if (!atSecECAL.isValid())
0826       continue;
0827     float secPhi = atSecECAL.positionREP().Phi();
0828
0829     unsigned int sbrem = 0;
0830     for (unsigned ipbrem = 0; ipbrem < primPFBrem.size(); ipbrem++) {
0831       if (primPFBrem[ipbrem].indTrajPoint() == 99)
0832         continue;
0833       const reco::PFTrajectoryPoint& atPrimECAL =
0834           primPFBrem[ipbrem].extrapolatedPoint(reco::PFTrajectoryPoint::ECALEntrance);
0835       if (!atPrimECAL.isValid())
0836         continue;
0837       sbrem++;
0838       if (sbrem <= 3) {
0839         float primPhi = atPrimECAL.positionREP().Phi();
0840
0841         float dphi = fabs(primPhi - secPhi);
0842         if (dphi > TMath::Pi())
0843           dphi -= TMath::TwoPi();
0844         if (fabs(dphi) < minDphi) {
0845           minDphi = fabs(dphi);
0846         }
0847       }
0848     }
0849
0850     cbrem++;
0851     if (cbrem == 3)
0852       break;
0853   }
0854   return minDphi;
0855 }
0856 bool PFElecTkProducer::isSameEgSC(const reco::ElectronSeed& nSeed,
0857                                   const reco::ElectronSeed& iSeed,
0858                                   bool& bothGsfEcalDriven,
0859                                   float& SCEnergy) {
0860   bool isSameSC = false;
0861
0862   if (nSeed.caloCluster().isNonnull() && iSeed.caloCluster().isNonnull()) {
0863     auto const* nscRef = static_cast<SuperCluster const*>(nSeed.caloCluster().get());
0864     auto const* iscRef = static_cast<SuperCluster const*>(iSeed.caloCluster().get());
0865
0866     if (nscRef && iscRef) {
0867       bothGsfEcalDriven = true;
0868       if (nscRef == iscRef) {
0869         isSameSC = true;
0870         // retrieve the supercluster energy
0871         SCEnergy = nscRef->energy();
0872       }
0873     }
0874   }
0875   return isSameSC;
0876 }
0877 bool PFElecTkProducer::isSharingEcalEnergyWithEgSC(const reco::GsfPFRecTrack& nGsfPFRecTrack,
0878                                                    const reco::GsfPFRecTrack& iGsfPFRecTrack,
0879                                                    const reco::ElectronSeed& nSeed,
0880                                                    const reco::ElectronSeed& iSeed,
0881                                                    const reco::PFClusterCollection& theEClus,
0882                                                    bool& bothGsfTrackerDriven,
0883                                                    bool& nEcalDriven,
0884                                                    bool& iEcalDriven,
0885                                                    float& nEnergy,
0886                                                    float& iEnergy) {
0887   bool isSharingEnergy = false;
0888
0889   //which is EcalDriven?
0890   bool oneEcalDriven = true;
0891   SuperCluster const* scRef = nullptr;
0892   GsfPFRecTrack gsfPfTrack;
0893
0894   if (nSeed.caloCluster().isNonnull()) {
0895     scRef = static_cast<SuperCluster const*>(nSeed.caloCluster().get());
0896     assert(scRef);
0897     nEnergy = scRef->energy();
0898     nEcalDriven = true;
0899     gsfPfTrack = iGsfPFRecTrack;
0900   } else if (iSeed.caloCluster().isNonnull()) {
0901     scRef = static_cast<SuperCluster const*>(iSeed.caloCluster().get());
0902     assert(scRef);
0903     iEnergy = scRef->energy();
0904     iEcalDriven = true;
0905     gsfPfTrack = nGsfPFRecTrack;
0906   } else {
0907     oneEcalDriven = false;
0908   }
0909
0910   if (oneEcalDriven) {
0911     //run a basic reconstruction for the particle flow
0912
0913     vector<PFCluster> vecPFClusters;
0914     vecPFClusters.clear();
0915
0916     for (PFClusterCollection::const_iterator clus = theEClus.begin(); clus != theEClus.end(); clus++) {
0917       PFCluster clust = *clus;
0918       clust.calculatePositionREP();
0919
0920       float deta = fabs(scRef->position().eta() - clust.position().eta());
0921       float dphi = fabs(scRef->position().phi() - clust.position().phi());
0922       if (dphi > TMath::Pi())
0923         dphi -= TMath::TwoPi();
0924
0925       // Angle preselection between the supercluster and pfclusters
0926       // this is needed just to save some cpu-time for this is hard-coded
0927       if (deta < 0.5 && fabs(dphi) < 1.0) {
0928         double distGsf = gsfPfTrack.extrapolatedPoint(reco::PFTrajectoryPoint::ECALShowerMax).isValid()
0929                              ? LinkByRecHit::testTrackAndClusterByRecHit(gsfPfTrack, clust)
0930                              : -1.;
0931         // check if it touch the GsfTrack
0932         if (distGsf > 0.) {
0933           if (nEcalDriven)
0934             iEnergy += clust.energy();
0935           else
0936             nEnergy += clust.energy();
0937           vecPFClusters.push_back(clust);
0938         }
0939         // check if it touch the Brem-tangents
0940         else {
0941           vector<PFBrem> primPFBrem = gsfPfTrack.PFRecBrem();
0942           for (unsigned ipbrem = 0; ipbrem < primPFBrem.size(); ipbrem++) {
0943             if (primPFBrem[ipbrem].indTrajPoint() == 99)
0944               continue;
0945             const reco::PFRecTrack& pfBremTrack = primPFBrem[ipbrem];
0946             double dist = pfBremTrack.extrapolatedPoint(reco::PFTrajectoryPoint::ECALShowerMax).isValid()
0947                               ? LinkByRecHit::testTrackAndClusterByRecHit(pfBremTrack, clust, true)
0948                               : -1.;
0949             if (dist > 0.) {
0950               if (nEcalDriven)
0951                 iEnergy += clust.energy();
0952               else
0953                 nEnergy += clust.energy();
0954               vecPFClusters.push_back(clust);
0955             }
0956           }
0957         }
0958       }  // END if angle preselection
0959     }    // PFClusters Loop
0960     if (!vecPFClusters.empty()) {
0961       for (unsigned int pf = 0; pf < vecPFClusters.size(); pf++) {
0962         bool isCommon = ClusterClusterMapping::overlap(vecPFClusters[pf], *scRef);
0963         if (isCommon) {
0964           isSharingEnergy = true;
0965         }
0966         break;
0967       }
0968     }
0969   } else {
0970     // both tracks are trackerDriven, try ECAL energy matching also in this case.
0971
0972     bothGsfTrackerDriven = true;
0973     vector<PFCluster> nPFCluster;
0974     vector<PFCluster> iPFCluster;
0975
0976     nPFCluster.clear();
0977     iPFCluster.clear();
0978
0979     for (PFClusterCollection::const_iterator clus = theEClus.begin(); clus != theEClus.end(); clus++) {
0980       PFCluster clust = *clus;
0981       clust.calculatePositionREP();
0982
0983       float ndeta = fabs(nGsfPFRecTrack.gsfTrackRef()->eta() - clust.position().eta());
0984       float ndphi = fabs(nGsfPFRecTrack.gsfTrackRef()->phi() - clust.position().phi());
0985       if (ndphi > TMath::Pi())
0986         ndphi -= TMath::TwoPi();
0987       // Apply loose preselection with the track
0988       // just to save cpu time, for this hard-coded
0989       if (ndeta < 0.5 && fabs(ndphi) < 1.0) {
0990         double distGsf = nGsfPFRecTrack.extrapolatedPoint(reco::PFTrajectoryPoint::ECALShowerMax).isValid()
0991                              ? LinkByRecHit::testTrackAndClusterByRecHit(nGsfPFRecTrack, clust)
0992                              : -1.;
0993         if (distGsf > 0.) {
0994           nPFCluster.push_back(clust);
0995           nEnergy += clust.energy();
0996         } else {
0997           const vector<PFBrem>& primPFBrem = nGsfPFRecTrack.PFRecBrem();
0998           for (unsigned ipbrem = 0; ipbrem < primPFBrem.size(); ipbrem++) {
0999             if (primPFBrem[ipbrem].indTrajPoint() == 99)
1000               continue;
1001             const reco::PFRecTrack& pfBremTrack = primPFBrem[ipbrem];
1002             double dist = pfBremTrack.extrapolatedPoint(reco::PFTrajectoryPoint::ECALShowerMax).isValid()
1003                               ? LinkByRecHit::testTrackAndClusterByRecHit(pfBremTrack, clust, true)
1004                               : -1.;
1005             if (dist > 0.) {
1006               nPFCluster.push_back(clust);
1007               nEnergy += clust.energy();
1008               break;
1009             }
1010           }
1011         }
1012       }
1013
1014       float ideta = fabs(iGsfPFRecTrack.gsfTrackRef()->eta() - clust.position().eta());
1015       float idphi = fabs(iGsfPFRecTrack.gsfTrackRef()->phi() - clust.position().phi());
1016       if (idphi > TMath::Pi())
1017         idphi -= TMath::TwoPi();
1018       // Apply loose preselection with the track
1019       // just to save cpu time, for this hard-coded
1020       if (ideta < 0.5 && fabs(idphi) < 1.0) {
1021         double dist = iGsfPFRecTrack.extrapolatedPoint(reco::PFTrajectoryPoint::ECALShowerMax).isValid()
1022                           ? LinkByRecHit::testTrackAndClusterByRecHit(iGsfPFRecTrack, clust)
1023                           : -1.;
1024         if (dist > 0.) {
1025           iPFCluster.push_back(clust);
1026           iEnergy += clust.energy();
1027         } else {
1028           vector<PFBrem> primPFBrem = iGsfPFRecTrack.PFRecBrem();
1029           for (unsigned ipbrem = 0; ipbrem < primPFBrem.size(); ipbrem++) {
1030             if (primPFBrem[ipbrem].indTrajPoint() == 99)
1031               continue;
1032             const reco::PFRecTrack& pfBremTrack = primPFBrem[ipbrem];
1033             double dist = LinkByRecHit::testTrackAndClusterByRecHit(pfBremTrack, clust, true);
1034             if (dist > 0.) {
1035               iPFCluster.push_back(clust);
1036               iEnergy += clust.energy();
1037               break;
1038             }
1039           }
1040         }
1041       }
1042     }
1043
1044     if (!nPFCluster.empty() && !iPFCluster.empty()) {
1045       for (unsigned int npf = 0; npf < nPFCluster.size(); npf++) {
1046         for (unsigned int ipf = 0; ipf < iPFCluster.size(); ipf++) {
1047           bool isCommon = ClusterClusterMapping::overlap(nPFCluster[npf], iPFCluster[ipf]);
1048           if (isCommon) {
1049             isSharingEnergy = true;
1050             break;
1051           }
1052         }
1053         if (isSharingEnergy)
1054           break;
1055       }
1056     }
1057   }
1058
1059   return isSharingEnergy;
1060 }
1061 bool PFElecTkProducer::isInnerMost(const reco::GsfTrackRef& nGsfTrack,
1062                                    const reco::GsfTrackRef& iGsfTrack,
1063                                    bool& sameLayer) {
1064   // copied by the class RecoEgamma/EgammaElectronAlgos/src/EgAmbiguityTools.cc
1065   // obsolete but the code is kept: now using lost hits method
1066
1067   const reco::HitPattern& gsfHitPattern1 = nGsfTrack->hitPattern();
1068   const reco::HitPattern& gsfHitPattern2 = iGsfTrack->hitPattern();
1069
1070   // retrieve first valid hit
1071   int gsfHitCounter1 = 0;
1072   for (auto const& hit : nGsfTrack->recHits()) {
1073     if (hit->isValid())
1074       break;
1075     gsfHitCounter1++;
1076   }
1077
1078   int gsfHitCounter2 = 0;
1079   for (auto const& hit : iGsfTrack->recHits()) {
1080     if (hit->isValid())
1081       break;
1082     gsfHitCounter2++;
1083   }
1084
1085   uint32_t gsfHit1 = gsfHitPattern1.getHitPattern(HitPattern::TRACK_HITS, gsfHitCounter1);
1086   uint32_t gsfHit2 = gsfHitPattern2.getHitPattern(HitPattern::TRACK_HITS, gsfHitCounter2);
1087
1088   if (gsfHitPattern1.getSubStructure(gsfHit1) != gsfHitPattern2.getSubStructure(gsfHit2)) {
1089     return (gsfHitPattern2.getSubStructure(gsfHit2) < gsfHitPattern1.getSubStructure(gsfHit1));
1090   } else if (gsfHitPattern1.getLayer(gsfHit1) != gsfHitPattern2.getLayer(gsfHit2)) {
1091     return (gsfHitPattern2.getLayer(gsfHit2) < gsfHitPattern1.getLayer(gsfHit1));
1092   } else {
1093     sameLayer = true;
1094     return false;
1095   }
1096 }
1097 bool PFElecTkProducer::isInnerMostWithLostHits(const reco::GsfTrackRef& nGsfTrack,
1098                                                const reco::GsfTrackRef& iGsfTrack,
1099                                                bool& sameLayer) {
1100   // define closest using the lost hits on the expectedhitsineer
1101   unsigned int nLostHits = nGsfTrack->missingInnerHits();
1102   unsigned int iLostHits = iGsfTrack->missingInnerHits();
1103
1104   if (nLostHits != iLostHits) {
1105     return (nLostHits > iLostHits);
1106   } else {
1107     sameLayer = true;
1108     return false;
1109   }
1110 }
1111
1112 // ------------ method called once each job just before starting event loop  ------------
1113 void PFElecTkProducer::beginRun(const edm::Run& run, const EventSetup& iSetup) {
1114   auto const& magneticField = &iSetup.getData(magFieldToken_);
1115   auto const& tracker = &iSetup.getData(tkerGeomToken_);
1116
1117   mtsTransform_ = MultiTrajectoryStateTransform(tracker, magneticField);
1118
1119   pfTransformer_ = std::make_unique<PFTrackTransformer>(math::XYZVector(magneticField->inTesla(GlobalPoint(0, 0, 0))));
1120
1121   TransientTrackBuilder thebuilder = iSetup.getData(transientTrackToken_);
1122
1123   convBremFinder_ = std::make_unique<ConvBremPFTrackFinder>(thebuilder,
1124                                                             mvaConvBremFinderIDBarrelLowPt_,
1125                                                             mvaConvBremFinderIDBarrelHighPt_,
1126                                                             mvaConvBremFinderIDEndcapsLowPt_,
1127                                                             mvaConvBremFinderIDEndcapsHighPt_);
1128 }
1129
1130 // ------------ method called once each job just after ending the event loop  ------------
1131 void PFElecTkProducer::endRun(const edm::Run& run, const EventSetup& iSetup) {
1132   pfTransformer_.reset();
1133   convBremFinder_.reset();
1134 }
1135
1136 //define this as a plug-in
1137 #include "FWCore/Framework/interface/MakerMacros.h"
1138 DEFINE_FWK_MODULE(PFElecTkProducer);