Back to home page

Project CMSSW displayed by LXR
 
 

    

File indexing completed on 2024-06-22 02:24:09 

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

This page was automatically generated by the 2.2.1 LXR engine.

The LXR team
Valid CSS 2.1! Valid HTML 4.01!