From 5c3f87339fb6c0af072da907f13e747cde001d0b Mon Sep 17 00:00:00 2001 From: Sam Harper Date: Thu, 31 Mar 2016 12:28:52 +0100 Subject: [PATCH 01/18] removing jetcoreregional tracks from tracker isolation sum --- .../interface/ElectronTkIsolation.h | 9 +++-- .../src/ElectronTkIsolation.cc | 35 ++++++++++++------- 2 files changed, 28 insertions(+), 16 deletions(-) diff --git a/RecoEgamma/EgammaIsolationAlgos/interface/ElectronTkIsolation.h b/RecoEgamma/EgammaIsolationAlgos/interface/ElectronTkIsolation.h index 3add1085d337b..915124322a5f9 100644 --- a/RecoEgamma/EgammaIsolationAlgos/interface/ElectronTkIsolation.h +++ b/RecoEgamma/EgammaIsolationAlgos/interface/ElectronTkIsolation.h @@ -45,7 +45,7 @@ class ElectronTkIsolation { drb_(drb), trackCollection_(trackCollection), beamPoint_(beamPoint) { - + setAlgosToReject(); setDzOption("vz"); } @@ -70,7 +70,7 @@ class ElectronTkIsolation { drb_(drb), trackCollection_(trackCollection), beamPoint_(beamPoint) { - + setAlgosToReject(); setDzOption("vz"); } @@ -106,7 +106,9 @@ class ElectronTkIsolation { std::pairgetIso(const reco::Track*) const ; private: - + + bool passAlgo(const reco::TrackBase& trk)const; + void setAlgosToReject(); double extRadius_ ; double intRadiusBarrel_ ; double intRadiusEndcap_ ; @@ -115,6 +117,7 @@ class ElectronTkIsolation { double ptLow_ ; double lip_ ; double drb_; + std::vector algosToReject_; //vector is sorted const reco::TrackCollection *trackCollection_ ; reco::TrackBase::Point beamPoint_; diff --git a/RecoEgamma/EgammaIsolationAlgos/src/ElectronTkIsolation.cc b/RecoEgamma/EgammaIsolationAlgos/src/ElectronTkIsolation.cc index 92c1eed3fead0..b957ddc69b1cb 100644 --- a/RecoEgamma/EgammaIsolationAlgos/src/ElectronTkIsolation.cc +++ b/RecoEgamma/EgammaIsolationAlgos/src/ElectronTkIsolation.cc @@ -45,6 +45,7 @@ ElectronTkIsolation::ElectronTkIsolation (double extRadius, trackCollection_(trackCollection), beamPoint_(beamPoint) { + setAlgosToReject(); setDzOption(dzOptionString); } @@ -86,19 +87,13 @@ std::pair ElectronTkIsolation::getIso(const reco::Track* tmpTrack) c if (fabs( (*itrTr).dxy(beamPoint_) ) > drb_ ) continue; double dr = ROOT::Math::VectorUtil::DeltaR(itrTr->momentum(),tmpElectronMomentumAtVtx) ; double deta = (*itrTr).eta() - tmpElectronEtaAtVertex; - if (fabs(tmpElectronEtaAtVertex) < 1.479) { - if ( fabs(dr) < extRadius_ && fabs(dr) >= intRadiusBarrel_ && fabs(deta) >= stripBarrel_) - { - ++counter ; - ptSum += this_pt; - } - } - else { - if ( fabs(dr) < extRadius_ && fabs(dr) >= intRadiusEndcap_ && fabs(deta) >= stripEndcap_) - { - ++counter ; - ptSum += this_pt; - } + bool isBarrel = std::abs(tmpElectronEtaAtVertex) < 1.479; + double intRadius = isBarrel ? intRadiusBarrel_ : intRadiusEndcap_; + double strip = isBarrel ? stripBarrel_ : stripEndcap_; + if(dr < extRadius_ && dr>=intRadius && std::abs(deta) >=strip && passAlgo(*itrTr)){ + + ++counter ; + ptSum += this_pt; } }//end loop over tracks @@ -122,3 +117,17 @@ double ElectronTkIsolation::getPtTracks (const reco::GsfElectron* electron) cons return getIso(electron).second ; } + +bool ElectronTkIsolation::passAlgo(const reco::TrackBase& trk)const +{ + int algo = trk.algo(); + bool rejAlgo=std::binary_search(algosToReject_.begin(),algosToReject_.end(),algo); + return rejAlgo==false; +} + + +void ElectronTkIsolation::setAlgosToReject() +{ + algosToReject_={reco::TrackBase::jetCoreRegionalStep}; + std::sort(algosToReject_.begin(),algosToReject_.end()); +} From b57291451df21d5ed8822d6e2b0ad17732ea788a Mon Sep 17 00:00:00 2001 From: Sam Harper Date: Sun, 5 Jun 2016 15:59:41 +0100 Subject: [PATCH 02/18] making EleTrk isolation have the option to set which algos to reject --- .../EgammaIsolationAlgos/interface/ElectronTkIsolation.h | 8 +++++--- .../EgammaIsolationAlgos/src/ElectronTkIsolation.cc | 7 ++++--- 2 files changed, 9 insertions(+), 6 deletions(-) diff --git a/RecoEgamma/EgammaIsolationAlgos/interface/ElectronTkIsolation.h b/RecoEgamma/EgammaIsolationAlgos/interface/ElectronTkIsolation.h index 915124322a5f9..350aafd678c10 100644 --- a/RecoEgamma/EgammaIsolationAlgos/interface/ElectronTkIsolation.h +++ b/RecoEgamma/EgammaIsolationAlgos/interface/ElectronTkIsolation.h @@ -45,7 +45,7 @@ class ElectronTkIsolation { drb_(drb), trackCollection_(trackCollection), beamPoint_(beamPoint) { - setAlgosToReject(); + setDefaultAlgosToReject(); setDzOption("vz"); } @@ -70,7 +70,7 @@ class ElectronTkIsolation { drb_(drb), trackCollection_(trackCollection), beamPoint_(beamPoint) { - setAlgosToReject(); + setDefaultAlgosToReject(); setDzOption("vz"); } @@ -104,11 +104,13 @@ class ElectronTkIsolation { double getPtTracks (const reco::GsfElectron*) const ; std::pairgetIso(const reco::GsfElectron*) const; std::pairgetIso(const reco::Track*) const ; + + void setAlgosToReject(std::vector algos); //intentional by value pass private: bool passAlgo(const reco::TrackBase& trk)const; - void setAlgosToReject(); + void setDefaultAlgosToReject(){} double extRadius_ ; double intRadiusBarrel_ ; double intRadiusEndcap_ ; diff --git a/RecoEgamma/EgammaIsolationAlgos/src/ElectronTkIsolation.cc b/RecoEgamma/EgammaIsolationAlgos/src/ElectronTkIsolation.cc index b957ddc69b1cb..f52f55e6b1680 100644 --- a/RecoEgamma/EgammaIsolationAlgos/src/ElectronTkIsolation.cc +++ b/RecoEgamma/EgammaIsolationAlgos/src/ElectronTkIsolation.cc @@ -45,7 +45,7 @@ ElectronTkIsolation::ElectronTkIsolation (double extRadius, trackCollection_(trackCollection), beamPoint_(beamPoint) { - setAlgosToReject(); + setDefaultAlgosToReject(); setDzOption(dzOptionString); } @@ -126,8 +126,9 @@ bool ElectronTkIsolation::passAlgo(const reco::TrackBase& trk)const } -void ElectronTkIsolation::setAlgosToReject() +void ElectronTkIsolation::setAlgosToReject(std::vector algos) { - algosToReject_={reco::TrackBase::jetCoreRegionalStep}; + algosToReject_=std::move(algos); std::sort(algosToReject_.begin(),algosToReject_.end()); } + From c858424218724780ec6688221cb96cbbef2a7382 Mon Sep 17 00:00:00 2001 From: Sam Harper Date: Mon, 6 Jun 2016 09:30:07 +0100 Subject: [PATCH 03/18] adding function to calculate number of saturated crystals in 5x5 --- .../interface/EcalClusterTools.h | 32 ++++++++++++++++++- 1 file changed, 31 insertions(+), 1 deletion(-) diff --git a/RecoEcal/EgammaCoreTools/interface/EcalClusterTools.h b/RecoEcal/EgammaCoreTools/interface/EcalClusterTools.h index 9d08200ce463b..5548c9f7fae1b 100644 --- a/RecoEcal/EgammaCoreTools/interface/EcalClusterTools.h +++ b/RecoEcal/EgammaCoreTools/interface/EcalClusterTools.h @@ -191,6 +191,9 @@ class EcalClusterToolsT { static std::vector roundnessBarrelSuperClusters( const reco::SuperCluster &superCluster ,const EcalRecHitCollection &recHits, int weightedPositionMethod = 0, float energyThreshold = 0.0); static std::vector roundnessBarrelSuperClustersUserExtended( const reco::SuperCluster &superCluster ,const EcalRecHitCollection &recHits, int ieta_delta=0, int iphi_delta=0, float energyRHThresh=0.00000, int weightedPositionMethod=0); static std::vector roundnessSelectedBarrelRecHits(const std::vector >&rhVector, int weightedPositionMethod = 0); + + //works out the number of staturated crystals in 5x5 + static int nrSaturatedCrysIn5x5(const DetId& id,const EcalRecHitCollection* recHits,const CaloTopology *topology); private: struct EcalClusterEnergyDeposition { @@ -247,7 +250,8 @@ class EcalClusterToolsT { static std::vector getSeedPosition(const std::vector >&RH_ptrs); static float getSumEnergy(const std::vector >&RH_ptrs_fracs); static float computeWeight(float eRH, float energyTotal, int weightedPositionMethod); - + + }; // implementation @@ -1626,6 +1630,31 @@ std::vector EcalClusterToolsT::roundnessSelectedBarrelRecHits( cons return shapes; } + + +template +int EcalClusterToolsT::nrSaturatedCrysIn5x5(const DetId& id,const EcalRecHitCollection* recHits,const CaloTopology *topology) +{ + int nrSat=0; + CaloNavigator cursor = CaloNavigator( id, topology->getSubdetectorTopology( id ) ); + + for ( int eastNr = -2; eastNr <= 2; ++eastNr ) { //east is eta in barrel + for ( int northNr = -2; northNr <= 2; ++northNr ) { //north is phi in barrel + cursor.home(); + cursor.offsetBy( eastNr, northNr); + DetId id = *cursor; + auto recHitIt = recHits->find(id); + if(recHitIt!=recHits->end() && + recHitIt->checkFlag(EcalRecHit::kSaturated)){ + nrSat++; + } + + } + } + return nrSat; +} + + //private functions useful for roundnessBarrelSuperClusters etc. //compute delta iphi between a seed and a particular recHit //iphi [1,360] @@ -1688,6 +1717,7 @@ float EcalClusterToolsT::getSumEnergy(const std::vector EcalClusterTools; namespace noZS { From 6c63e265aa6b208b5f8f6c7b21a4f28639bb317d Mon Sep 17 00:00:00 2001 From: Sam Harper Date: Mon, 6 Jun 2016 09:30:32 +0100 Subject: [PATCH 04/18] add value map producer for HEEP ID --- .../ElectronIdentification/BuildFile.xml | 1 + .../plugins/ElectronHEEPIDValueMapProducer.cc | 181 ++++++++++++++++++ 2 files changed, 182 insertions(+) create mode 100644 RecoEgamma/ElectronIdentification/plugins/ElectronHEEPIDValueMapProducer.cc diff --git a/RecoEgamma/ElectronIdentification/BuildFile.xml b/RecoEgamma/ElectronIdentification/BuildFile.xml index 461818d4075fd..f7589cf2db2d7 100644 --- a/RecoEgamma/ElectronIdentification/BuildFile.xml +++ b/RecoEgamma/ElectronIdentification/BuildFile.xml @@ -7,6 +7,7 @@ + diff --git a/RecoEgamma/ElectronIdentification/plugins/ElectronHEEPIDValueMapProducer.cc b/RecoEgamma/ElectronIdentification/plugins/ElectronHEEPIDValueMapProducer.cc new file mode 100644 index 0000000000000..b80626bd2e541 --- /dev/null +++ b/RecoEgamma/ElectronIdentification/plugins/ElectronHEEPIDValueMapProducer.cc @@ -0,0 +1,181 @@ +#include "FWCore/Framework/interface/Frameworkfwd.h" +#include "FWCore/Framework/interface/stream/EDProducer.h" + +#include "FWCore/Framework/interface/Event.h" +#include "FWCore/Framework/interface/MakerMacros.h" + +#include "FWCore/ParameterSet/interface/ParameterSet.h" + +#include "DataFormats/Common/interface/ValueMap.h" +#include "DataFormats/Common/interface/View.h" + +#include "DataFormats/EgammaCandidates/interface/GsfElectron.h" +#include "DataFormats/EgammaCandidates/interface/GsfElectronFwd.h" +#include "DataFormats/EcalDetId/interface/EcalSubdetector.h" + +#include "Geometry/CaloTopology/interface/CaloTopology.h" +#include "Geometry/CaloEventSetup/interface/CaloTopologyRecord.h" + +#include "RecoEcal/EgammaCoreTools/interface/EcalClusterTools.h" +#include "RecoEgamma/EgammaIsolationAlgos/interface/ElectronTkIsolation.h" + +#include +#include + +//Heavily inspired from ElectronIDValueMapProducer + + +class ElectronHEEPIDValueMapProducer : public edm::stream::EDProducer<> { + + public: + + explicit ElectronHEEPIDValueMapProducer(const edm::ParameterSet&); + ~ElectronHEEPIDValueMapProducer(); + + static void fillDescriptions(edm::ConfigurationDescriptions& descriptions); + + private: + + virtual void produce(edm::Event&, const edm::EventSetup&) override; + + template + static void writeValueMap(edm::Event &iEvent, + const edm::Handle > & handle, + const std::vector & values, + const std::string& label); + + static int nrSaturatedCrysIn5x5(const reco::GsfElectron& ele, + edm::Handle& ebHits, + edm::Handle& eeHits, + edm::ESHandle& caloTopo); + + template void setToken(edm::EDGetTokenT& token,edm::InputTag tag){token=consumes(tag);} + template void setToken(edm::EDGetTokenT& token,const edm::ParameterSet& iPara,const std::string& tag){token=consumes(iPara.getParameter(tag));} + template edm::Handle getHandle(const edm::Event& iEvent,const edm::EDGetTokenT& token){ + edm::Handle handle; + iEvent.getByToken(token,handle); + return handle; + } + + + edm::EDGetTokenT ebRecHitToken_; + edm::EDGetTokenT eeRecHitToken_; + edm::EDGetTokenT > eleToken_; + edm::EDGetTokenT beamSpotToken_; + edm::EDGetTokenT trkToken_; + + struct TrkIsoParam { + double extRadius; + double intRadiusBarrel; + double intRadiusEndcap; + double stripBarrel; + double stripEndcap; + double ptLow ; + double lip ; + double drb; + TrkIsoParam(const edm::ParameterSet& iPara); + }; + const TrkIsoParam trkIsoParam_; + static const std::string eleTrkPtIsoNoJetCoreLabel_; + static const std::string eleNrSaturateIn5x5Label_; +}; + +const std::string ElectronHEEPIDValueMapProducer::eleTrkPtIsoNoJetCoreLabel_="eleTrkPtIsoNoJetCore"; +const std::string ElectronHEEPIDValueMapProducer::eleNrSaturateIn5x5Label_="eleNrSaturateIn5x5"; + +ElectronHEEPIDValueMapProducer::TrkIsoParam::TrkIsoParam(const edm::ParameterSet& iParam): + extRadius(iParam.getParameter("extRadius")), + intRadiusBarrel(iParam.getParameter("intRadiusBarrel")), + intRadiusEndcap(iParam.getParameter("intRadiusEndcap")), + stripBarrel(iParam.getParameter("stripBarrel")), + stripEndcap(iParam.getParameter("stripEndcap")), + ptLow(iParam.getParameter("ptLow")), + lip(iParam.getParameter("lip")), + drb(iParam.getParameter("drb")) +{ + +} + + +ElectronHEEPIDValueMapProducer::ElectronHEEPIDValueMapProducer(const edm::ParameterSet& iConfig): + trkIsoParam_(iConfig.getParameter("trkIsoConfig")) +{ + setToken(ebRecHitToken_,iConfig,"ebRecHits"); + setToken(eeRecHitToken_,iConfig,"eeRecHits"); + setToken(eleToken_,iConfig,"eles"); + setToken(trkToken_,iConfig,"tracks"); + setToken(beamSpotToken_,iConfig,"beamSpot"); + + produces >(eleTrkPtIsoNoJetCoreLabel_); + produces >(eleNrSaturateIn5x5Label_); +} + +ElectronHEEPIDValueMapProducer::~ElectronHEEPIDValueMapProducer() +{ + +} + +void ElectronHEEPIDValueMapProducer::produce(edm::Event& iEvent, const edm::EventSetup& iSetup) +{ + auto eleHandle = getHandle(iEvent,eleToken_); + auto trkHandle = getHandle(iEvent,trkToken_); + auto ebRecHitHandle = getHandle(iEvent,ebRecHitToken_); + auto eeRecHitHandle = getHandle(iEvent,eeRecHitToken_); + auto beamSpotHandle = getHandle(iEvent,beamSpotToken_); + + edm::ESHandle caloTopoHandle; + iSetup.get().get(caloTopoHandle); + + ElectronTkIsolation isolCorr(trkIsoParam_.extRadius,trkIsoParam_.intRadiusBarrel,trkIsoParam_.intRadiusEndcap, + trkIsoParam_.stripBarrel,trkIsoParam_.stripEndcap,trkIsoParam_.ptLow, + trkIsoParam_.lip,trkIsoParam_.drb, + trkHandle.product(),beamSpotHandle->position()); + isolCorr.setAlgosToReject({reco::TrackBase::jetCoreRegionalStep}); + + std::vector eleTrkPtIsoNoJetCore; + std::vector eleNrSaturateIn5x5; + + for(size_t eleNr=0;eleNrsize();eleNr++){ + auto elePtr = eleHandle->ptrAt(eleNr); + eleTrkPtIsoNoJetCore.push_back(isolCorr.getPtTracks(&*elePtr)); + eleNrSaturateIn5x5.push_back(nrSaturatedCrysIn5x5(*elePtr,ebRecHitHandle,eeRecHitHandle,caloTopoHandle)); + } + + writeValueMap(iEvent,eleHandle,eleTrkPtIsoNoJetCore,eleTrkPtIsoNoJetCoreLabel_); + writeValueMap(iEvent,eleHandle,eleNrSaturateIn5x5,eleNrSaturateIn5x5Label_); + +} + +int ElectronHEEPIDValueMapProducer::nrSaturatedCrysIn5x5(const reco::GsfElectron& ele, + edm::Handle& ebHits, + edm::Handle& eeHits, + edm::ESHandle& caloTopo) +{ + DetId id = ele.superCluster()->seed()->seed(); + auto recHits = id.subdetId()==EcalBarrel ? ebHits.product() : eeHits.product(); + return noZS::EcalClusterTools::nrSaturatedCrysIn5x5(id,recHits,caloTopo.product()); + +} + +template +void ElectronHEEPIDValueMapProducer::writeValueMap(edm::Event &iEvent, + const edm::Handle > & handle, + const std::vector & values, + const std::string& label) +{ + std::unique_ptr > valMap(new edm::ValueMap()); + typename edm::ValueMap::Filler filler(*valMap); + filler.insert(handle, values.begin(), values.end()); + filler.fill(); + iEvent.put(std::move(valMap),label); +} + +void ElectronHEEPIDValueMapProducer::fillDescriptions(edm::ConfigurationDescriptions& descriptions) { + //The following says we do not know what parameters are allowed so do no validation + // Please change this to state exactly what you do use, even if it is no parameters + edm::ParameterSetDescription desc; + desc.setUnknown(); + descriptions.addDefault(desc); +} + +DEFINE_FWK_MODULE(ElectronHEEPIDValueMapProducer); From c6bf88063e9adcc901a09f1802847ca664b4e197 Mon Sep 17 00:00:00 2001 From: Sam Harper Date: Fri, 15 Apr 2016 16:05:21 +0100 Subject: [PATCH 05/18] adding in HEEP V6.1 patch --- .../plugins/cuts/GsfEleTrkPtIsoRhoCut.cc | 80 +++++++++++++ .../Identification/heepElectronID_tools.py | 109 ++++++++++++++++++ 2 files changed, 189 insertions(+) create mode 100644 RecoEgamma/ElectronIdentification/plugins/cuts/GsfEleTrkPtIsoRhoCut.cc diff --git a/RecoEgamma/ElectronIdentification/plugins/cuts/GsfEleTrkPtIsoRhoCut.cc b/RecoEgamma/ElectronIdentification/plugins/cuts/GsfEleTrkPtIsoRhoCut.cc new file mode 100644 index 0000000000000..f99a133d0c86d --- /dev/null +++ b/RecoEgamma/ElectronIdentification/plugins/cuts/GsfEleTrkPtIsoRhoCut.cc @@ -0,0 +1,80 @@ +#include "PhysicsTools/SelectorUtils/interface/CutApplicatorWithEventContentBase.h" +#include "DataFormats/EgammaCandidates/interface/GsfElectron.h" +#include "DataFormats/EgammaCandidates/interface/ConversionFwd.h" +#include "DataFormats/EgammaCandidates/interface/Conversion.h" +#include "RecoEgamma/EgammaTools/interface/ConversionTools.h" + +#include "RecoEgamma/ElectronIdentification/interface/EBEECutValues.h" + +class GsfEleTrkPtIsoRhoCut : public CutApplicatorWithEventContentBase { +public: + GsfEleTrkPtIsoRhoCut(const edm::ParameterSet& c); + + result_type operator()(const reco::GsfElectronPtr&) const override final; + + void setConsumes(edm::ConsumesCollector&) override final; + void getEventContent(const edm::EventBase&) override final; + + double value(const reco::CandidatePtr& cand) const override final; + + CandidateType candidateType() const override final { + return ELECTRON; + } + +private: + + EBEECutValues slopeTerm_; + EBEECutValues slopeStart_; + EBEECutValues constTerm_; + EBEECutValues rhoEtStart_; + EBEECutValues rhoEA_; + + + edm::Handle rhoHandle_; + +}; + +DEFINE_EDM_PLUGIN(CutApplicatorFactory, + GsfEleTrkPtIsoRhoCut, + "GsfEleTrkPtIsoRhoCut"); + +GsfEleTrkPtIsoRhoCut::GsfEleTrkPtIsoRhoCut(const edm::ParameterSet& params) : + CutApplicatorWithEventContentBase(params), + slopeTerm_(params,"slopeTerm"), + slopeStart_(params,"slopeStart"), + constTerm_(params,"constTerm"), + rhoEtStart_(params,"rhoEtStart"), + rhoEA_(params,"rhoEA") +{ + edm::InputTag rhoTag = params.getParameter("rho"); + contentTags_.emplace("rho",rhoTag); +} + +void GsfEleTrkPtIsoRhoCut::setConsumes(edm::ConsumesCollector& cc) { + auto rho = cc.consumes(contentTags_["rho"]); + contentTokens_.emplace("rho",rho); +} + +void GsfEleTrkPtIsoRhoCut::getEventContent(const edm::EventBase& ev) { + ev.getByLabel(contentTags_["rho"],rhoHandle_); +} + +CutApplicatorBase::result_type +GsfEleTrkPtIsoRhoCut:: +operator()(const reco::GsfElectronPtr& cand) const{ + const double rho = (*rhoHandle_); + const float isolTrkPt = cand->dr03TkSumPt(); + + const float et = cand->et(); + const float cutValue = et > slopeStart_(cand) ? slopeTerm_(cand)*(et-slopeStart_(cand)) + constTerm_(cand) : constTerm_(cand); + + const float rhoCutValue = et >= rhoEtStart_(cand) ? rhoEA_(cand)*rho : 0.; + + return isolTrkPt < cutValue+rhoCutValue; +} + +double GsfEleTrkPtIsoRhoCut:: +value(const reco::CandidatePtr& cand) const { + reco::GsfElectronPtr ele(cand); + return ele->dr03TkSumPt(); +} diff --git a/RecoEgamma/ElectronIdentification/python/Identification/heepElectronID_tools.py b/RecoEgamma/ElectronIdentification/python/Identification/heepElectronID_tools.py index afb49da39508b..8fabf19234354 100644 --- a/RecoEgamma/ElectronIdentification/python/Identification/heepElectronID_tools.py +++ b/RecoEgamma/ElectronIdentification/python/Identification/heepElectronID_tools.py @@ -68,6 +68,69 @@ def __init__(self, self.dxyCut = dxyCut self.maxMissingHitsCut = maxMissingHitsCut self.ecalDrivenCut = ecalDrivenCut + +class HEEP_WorkingPoint_V2: + """ + This is a container class to hold numerical cut values for either + the barrel or endcap set of cuts + """ + def __init__(self, + idName, + dEtaInSeedCut, + dPhiInCut, + full5x5SigmaIEtaIEtaCut, + # Two constants for the GsfEleFull5x5E2x5OverE5x5Cut + minE1x5OverE5x5Cut, + minE2x5OverE5x5Cut, + # Three constants for the GsfEleHadronicOverEMLinearCut: + # cut = constTerm if value < slopeStart + # cut = slopeTerm * (value - slopeStart) + constTerm if value >= slopeStart + hOverESlopeTerm, + hOverESlopeStart, + hOverEConstTerm, + # Three constants for the GsfEleTrkPtIsoCut: + # cut = constTerm if value < slopeStart + # cut = slopeTerm * (value - slopeStart) + constTerm if value >= slopeStart + trkIsoSlopeTerm, + trkIsoSlopeStart, + trkIsoConstTerm, + trkIsoRhoCorrStart, + trkIsoEffArea, + # Three constants for the GsfEleEmHadD1IsoRhoCut: + # cut = constTerm if value < slopeStart + # cut = slopeTerm * (value - slopeStart) + constTerm if value >= slopeStart + # Also for the same cut, the effective area for the rho correction of the isolation + ehIsoSlopeTerm, + ehIsoSlopeStart, + ehIsoConstTerm, + effAreaForEHIso, + # other cuts: + dxyCut, + maxMissingHitsCut, + ecalDrivenCut + ): + # assign values taken from all the arguments above + self.idName = idName + self.dEtaInSeedCut = dEtaInSeedCut + self.dPhiInCut = dPhiInCut + self.full5x5SigmaIEtaIEtaCut = full5x5SigmaIEtaIEtaCut + self.minE1x5OverE5x5Cut = minE1x5OverE5x5Cut + self.minE2x5OverE5x5Cut = minE2x5OverE5x5Cut + self.hOverESlopeTerm = hOverESlopeTerm + self.hOverESlopeStart = hOverESlopeStart + self.hOverEConstTerm = hOverEConstTerm + self.trkIsoSlopeTerm = trkIsoSlopeTerm + self.trkIsoSlopeStart = trkIsoSlopeStart + self.trkIsoConstTerm = trkIsoConstTerm + self.trkIsoRhoCorrStart = trkIsoRhoCorrStart + self.trkIsoEffArea = trkIsoEffArea + self.ehIsoSlopeTerm = ehIsoSlopeTerm + self.ehIsoSlopeStart = ehIsoSlopeStart + self.ehIsoConstTerm = ehIsoConstTerm + self.effAreaForEHIso = effAreaForEHIso + self.dxyCut = dxyCut + self.maxMissingHitsCut = maxMissingHitsCut + self.ecalDrivenCut = ecalDrivenCut # ============================================================== # Define individual cut configurations used by complete cut sets # ============================================================== @@ -176,6 +239,27 @@ def psetGsfEleTrkPtIsoCut(wpEB, wpEE): needsAdditionalProducts = cms.bool(False), isIgnored = cms.bool(False) ) +# Configure the cut on the tracker isolation with a rho correction (hack for 76X) +def psetGsfEleTrkPtIsoRhoCut(wpEB, wpEE): + return cms.PSet( + cutName = cms.string('GsfEleTrkPtIsoRhoCut'), + # Three constants for the GsfEleTrkPtIsoCut + # cut = constTerm if value < slopeStart + # cut = slopeTerm * (value - slopeStart) + constTerm if value >= slopeStart + slopeTermEB = cms.double( wpEB.trkIsoSlopeTerm ), + slopeTermEE = cms.double( wpEE.trkIsoSlopeTerm ), + slopeStartEB = cms.double( wpEB.trkIsoSlopeStart ), + slopeStartEE = cms.double( wpEE.trkIsoSlopeStart ), + constTermEB = cms.double( wpEB.trkIsoConstTerm ), + constTermEE = cms.double( wpEE.trkIsoConstTerm ), + rhoEtStartEB = cms.double( wpEB.trkIsoRhoCorrStart), + rhoEtStartEE = cms.double( wpEE.trkIsoRhoCorrStart), + rhoEAEB = cms.double( wpEB.trkIsoEffArea), + rhoEAEE = cms.double( wpEE.trkIsoEffArea), + rho = cms.InputTag("fixedGridRhoFastjetAll"), + needsAdditionalProducts = cms.bool(True), + isIgnored = cms.bool(False) + ) # Configure the cut on the EM + Had_depth_1 isolation with rho correction def psetGsfEleEmHadD1IsoRhoCut(wpEB, wpEE): @@ -279,3 +363,28 @@ def configureHEEPElectronID_V60(wpEB, wpEE): ) ) return parameterSet + +def configureHEEPElectronID_V61(wpEB, wpEE): + """ + This function configures the full cms.PSet for a VID ID and returns it. + The inputs: two objects of the type HEEP_WorkingPoint_V2, one + containing the cuts for the Barrel (EB) and the other one for the Endcap (EE). + """ + parameterSet = cms.PSet( + idName = cms.string("heepElectronID-HEEPV61"), + cutFlow = cms.VPSet( + psetMinPtCut(), #0 + psetGsfEleSCEtaMultiRangeCut(), #1 + psetGsfEleDEtaInSeedCut(wpEB,wpEE), #2 + psetGsfEleDPhiInCut(wpEB,wpEE), #3 + psetGsfEleFull5x5SigmaIEtaIEtaCut(wpEB,wpEE), #4 + psetGsfEleFull5x5E2x5OverE5x5Cut(wpEB,wpEE), #5 + psetGsfEleHadronicOverEMLinearCut(wpEB,wpEE), #6 + psetGsfEleTrkPtIsoRhoCut(wpEB,wpEE), #7 + psetGsfEleEmHadD1IsoRhoCut(wpEB,wpEE), #8 + psetGsfEleDxyCut(wpEB,wpEE), #9 + psetGsfEleMissingHitsCut(wpEB,wpEE), #10, + psetGsfEleEcalDrivenCut(wpEB,wpEE) #11 + ) + ) + return parameterSet From 644ffbdcf7dd17131bb730a925dbd54b463ed222 Mon Sep 17 00:00:00 2001 From: Sam Harper Date: Mon, 6 Jun 2016 12:11:16 +0100 Subject: [PATCH 06/18] changing track isol config name --- .../plugins/ElectronHEEPIDValueMapProducer.cc | 12 ++++++------ 1 file changed, 6 insertions(+), 6 deletions(-) diff --git a/RecoEgamma/ElectronIdentification/plugins/ElectronHEEPIDValueMapProducer.cc b/RecoEgamma/ElectronIdentification/plugins/ElectronHEEPIDValueMapProducer.cc index b80626bd2e541..a7028863dced3 100644 --- a/RecoEgamma/ElectronIdentification/plugins/ElectronHEEPIDValueMapProducer.cc +++ b/RecoEgamma/ElectronIdentification/plugins/ElectronHEEPIDValueMapProducer.cc @@ -70,8 +70,8 @@ class ElectronHEEPIDValueMapProducer : public edm::stream::EDProducer<> { double intRadiusEndcap; double stripBarrel; double stripEndcap; - double ptLow ; - double lip ; + double ptMin; + double maxVtxDist; double drb; TrkIsoParam(const edm::ParameterSet& iPara); }; @@ -89,8 +89,8 @@ ElectronHEEPIDValueMapProducer::TrkIsoParam::TrkIsoParam(const edm::ParameterSet intRadiusEndcap(iParam.getParameter("intRadiusEndcap")), stripBarrel(iParam.getParameter("stripBarrel")), stripEndcap(iParam.getParameter("stripEndcap")), - ptLow(iParam.getParameter("ptLow")), - lip(iParam.getParameter("lip")), + ptMin(iParam.getParameter("ptMin")), + maxVtxDist(iParam.getParameter("maxVtxDist")), drb(iParam.getParameter("drb")) { @@ -127,8 +127,8 @@ void ElectronHEEPIDValueMapProducer::produce(edm::Event& iEvent, const edm::Even iSetup.get().get(caloTopoHandle); ElectronTkIsolation isolCorr(trkIsoParam_.extRadius,trkIsoParam_.intRadiusBarrel,trkIsoParam_.intRadiusEndcap, - trkIsoParam_.stripBarrel,trkIsoParam_.stripEndcap,trkIsoParam_.ptLow, - trkIsoParam_.lip,trkIsoParam_.drb, + trkIsoParam_.stripBarrel,trkIsoParam_.stripEndcap,trkIsoParam_.ptMin, + trkIsoParam_.maxVtxDist,trkIsoParam_.drb, trkHandle.product(),beamSpotHandle->position()); isolCorr.setAlgosToReject({reco::TrackBase::jetCoreRegionalStep}); From 4e208401235dc79b66bafb1fb750a257c9353d38 Mon Sep 17 00:00:00 2001 From: Sam Harper Date: Mon, 6 Jun 2016 12:12:05 +0100 Subject: [PATCH 07/18] adding 1st version of HEEP V6.0 in 80X --- .../GsfEleFull5x5E2x5OverE5x5WithSatCut.cc | 69 ++++++++++++ .../GsfEleFull5x5SigmaIEtaIEtaWithSatCut.cc | 63 +++++++++++ .../plugins/cuts/GsfEleValueMapIsoRhoCut.cc | 94 ++++++++++++++++ .../heepElectronID_HEEPV60_80XAOD_cff.py | 102 ++++++++++++++++++ .../Identification/heepElectronID_tools.py | 76 ++++++++++++- 5 files changed, 403 insertions(+), 1 deletion(-) create mode 100644 RecoEgamma/ElectronIdentification/plugins/cuts/GsfEleFull5x5E2x5OverE5x5WithSatCut.cc create mode 100644 RecoEgamma/ElectronIdentification/plugins/cuts/GsfEleFull5x5SigmaIEtaIEtaWithSatCut.cc create mode 100644 RecoEgamma/ElectronIdentification/plugins/cuts/GsfEleValueMapIsoRhoCut.cc create mode 100644 RecoEgamma/ElectronIdentification/python/Identification/heepElectronID_HEEPV60_80XAOD_cff.py diff --git a/RecoEgamma/ElectronIdentification/plugins/cuts/GsfEleFull5x5E2x5OverE5x5WithSatCut.cc b/RecoEgamma/ElectronIdentification/plugins/cuts/GsfEleFull5x5E2x5OverE5x5WithSatCut.cc new file mode 100644 index 0000000000000..2ff6d17f4c258 --- /dev/null +++ b/RecoEgamma/ElectronIdentification/plugins/cuts/GsfEleFull5x5E2x5OverE5x5WithSatCut.cc @@ -0,0 +1,69 @@ +#include "PhysicsTools/SelectorUtils/interface/CutApplicatorWithEventContentBase.h" +#include "DataFormats/EgammaCandidates/interface/GsfElectron.h" +#include "DataFormats/Common/interface/ValueMap.h" +#include "RecoEgamma/ElectronIdentification/interface/EBEECutValues.h" + +class GsfEleFull5x5E2x5OverE5x5WithSatCut : public CutApplicatorWithEventContentBase { +public: + GsfEleFull5x5E2x5OverE5x5WithSatCut(const edm::ParameterSet& c); + + result_type operator()(const reco::GsfElectronPtr&) const override final; + + void setConsumes(edm::ConsumesCollector&) override final; + void getEventContent(const edm::EventBase&) override final; + + double value(const reco::CandidatePtr& cand) const override final; + + CandidateType candidateType() const override final { + return ELECTRON; + } + +private: + EBEECutValues minE1x5OverE5x5Cut_; + EBEECutValues minE2x5OverE5x5Cut_; + EBEECutValuesInt maxNrSatCrysIn5x5Cut_; + edm::Handle > nrSatCrysValueMap_; + +}; + +DEFINE_EDM_PLUGIN(CutApplicatorFactory, + GsfEleFull5x5E2x5OverE5x5WithSatCut, + "GsfEleFull5x5E2x5OverE5x5WithSatCut"); + +GsfEleFull5x5E2x5OverE5x5WithSatCut::GsfEleFull5x5E2x5OverE5x5WithSatCut(const edm::ParameterSet& params) : + CutApplicatorWithEventContentBase(params), + minE1x5OverE5x5Cut_(params,"minE1x5OverE5x5"), + minE2x5OverE5x5Cut_(params,"minE2x5OverE5x5"), + maxNrSatCrysIn5x5Cut_(params,"maxNrSatCrysIn5x5"){ + contentTags_.emplace("nrSatCrysValueMap",params.getParameter("nrSatCrysValueMap")); +} + +void GsfEleFull5x5E2x5OverE5x5WithSatCut::setConsumes(edm::ConsumesCollector& cc) { + contentTokens_.emplace("nrSatCrysValueMap",cc.consumes >(contentTags_["nrSatCrysValueMap"])); +} + +void GsfEleFull5x5E2x5OverE5x5WithSatCut::getEventContent(const edm::EventBase& ev) { + ev.getByLabel(contentTags_["nrSatCrysValueMap"],nrSatCrysValueMap_); + +} + + +CutApplicatorBase::result_type +GsfEleFull5x5E2x5OverE5x5WithSatCut:: +operator()(const reco::GsfElectronPtr& cand) const{ + + if((*nrSatCrysValueMap_)[cand]>maxNrSatCrysIn5x5Cut_(cand)) return false; + + const double e5x5 = cand->full5x5_e5x5(); + const double e2x5OverE5x5 = e5x5!=0 ? cand->full5x5_e2x5Max()/e5x5 : 0; + const double e1x5OverE5x5 = e5x5!=0 ? cand->full5x5_e1x5()/e5x5 : 0; + + return e1x5OverE5x5 > minE1x5OverE5x5Cut_(cand) || e2x5OverE5x5 > minE2x5OverE5x5Cut_(cand); + +} + +double GsfEleFull5x5E2x5OverE5x5WithSatCut:: +value(const reco::CandidatePtr& cand) const { + reco::GsfElectronPtr ele(cand); + return ele->full5x5_e2x5Max()/ele->full5x5_e1x5(); +} diff --git a/RecoEgamma/ElectronIdentification/plugins/cuts/GsfEleFull5x5SigmaIEtaIEtaWithSatCut.cc b/RecoEgamma/ElectronIdentification/plugins/cuts/GsfEleFull5x5SigmaIEtaIEtaWithSatCut.cc new file mode 100644 index 0000000000000..78d197966e905 --- /dev/null +++ b/RecoEgamma/ElectronIdentification/plugins/cuts/GsfEleFull5x5SigmaIEtaIEtaWithSatCut.cc @@ -0,0 +1,63 @@ +#include "PhysicsTools/SelectorUtils/interface/CutApplicatorWithEventContentBase.h" +#include "DataFormats/EgammaCandidates/interface/GsfElectron.h" +#include "DataFormats/Common/interface/ValueMap.h" +#include "RecoEgamma/ElectronIdentification/interface/EBEECutValues.h" + +class GsfEleFull5x5SigmaIEtaIEtaWithSatCut : public CutApplicatorWithEventContentBase { +public: + GsfEleFull5x5SigmaIEtaIEtaWithSatCut(const edm::ParameterSet& c); + + result_type operator()(const reco::GsfElectronPtr&) const override final; + + void setConsumes(edm::ConsumesCollector&) override final; + void getEventContent(const edm::EventBase&) override final; + + double value(const reco::CandidatePtr& cand) const override final; + + CandidateType candidateType() const override final { + return ELECTRON; + } + +private: + EBEECutValues maxSigmaIEtaIEtaCut_; + EBEECutValuesInt maxNrSatCrysIn5x5Cut_; + edm::Handle > nrSatCrysValueMap_; + +}; + +DEFINE_EDM_PLUGIN(CutApplicatorFactory, + GsfEleFull5x5SigmaIEtaIEtaWithSatCut, + "GsfEleFull5x5SigmaIEtaIEtaWithSatCut"); + +GsfEleFull5x5SigmaIEtaIEtaWithSatCut::GsfEleFull5x5SigmaIEtaIEtaWithSatCut(const edm::ParameterSet& params) : + CutApplicatorWithEventContentBase(params), + maxSigmaIEtaIEtaCut_(params,"maxSigmaIEtaIEta"), + + maxNrSatCrysIn5x5Cut_(params,"maxNrSatCrysIn5x5"){ + contentTags_.emplace("nrSatCrysValueMap",params.getParameter("nrSatCrysValueMap")); +} + +void GsfEleFull5x5SigmaIEtaIEtaWithSatCut::setConsumes(edm::ConsumesCollector& cc) { + contentTokens_.emplace("nrSatCrysValueMap",cc.consumes >(contentTags_["nrSatCrysValueMap"])); +} + +void GsfEleFull5x5SigmaIEtaIEtaWithSatCut::getEventContent(const edm::EventBase& ev) { + ev.getByLabel(contentTags_["nrSatCrysValueMap"],nrSatCrysValueMap_); + +} + + +CutApplicatorBase::result_type +GsfEleFull5x5SigmaIEtaIEtaWithSatCut:: +operator()(const reco::GsfElectronPtr& cand) const{ + + if((*nrSatCrysValueMap_)[cand]>maxNrSatCrysIn5x5Cut_(cand)) return false; + else return cand->full5x5_sigmaIetaIeta() < maxSigmaIEtaIEtaCut_(cand); + +} + +double GsfEleFull5x5SigmaIEtaIEtaWithSatCut:: +value(const reco::CandidatePtr& cand) const { + reco::GsfElectronPtr ele(cand); + return ele->full5x5_sigmaIetaIeta(); +} diff --git a/RecoEgamma/ElectronIdentification/plugins/cuts/GsfEleValueMapIsoRhoCut.cc b/RecoEgamma/ElectronIdentification/plugins/cuts/GsfEleValueMapIsoRhoCut.cc new file mode 100644 index 0000000000000..6d9f23a38f61c --- /dev/null +++ b/RecoEgamma/ElectronIdentification/plugins/cuts/GsfEleValueMapIsoRhoCut.cc @@ -0,0 +1,94 @@ +#include "PhysicsTools/SelectorUtils/interface/CutApplicatorWithEventContentBase.h" +#include "DataFormats/EgammaCandidates/interface/GsfElectron.h" +#include "DataFormats/Common/interface/ValueMap.h" +#include "RecoEgamma/ElectronIdentification/interface/EBEECutValues.h" + +class GsfEleValueMapIsoRhoCut : public CutApplicatorWithEventContentBase { +public: + GsfEleValueMapIsoRhoCut(const edm::ParameterSet& c); + + result_type operator()(const reco::GsfElectronPtr&) const override final; + + void setConsumes(edm::ConsumesCollector&) override final; + void getEventContent(const edm::EventBase&) override final; + + double value(const reco::CandidatePtr& cand) const override final; + + CandidateType candidateType() const override final { + return ELECTRON; + } + +private: + float getRhoCorr(const reco::GsfElectronPtr& cand)const; + + EBEECutValues slopeTerm_; + EBEECutValues slopeStart_; + EBEECutValues constTerm_; + EBEECutValues rhoEtStart_; + EBEECutValues rhoEA_; + + bool useRho_; + + edm::Handle rhoHandle_; + edm::Handle > valueHandle_; +}; + +DEFINE_EDM_PLUGIN(CutApplicatorFactory, + GsfEleValueMapIsoRhoCut, + "GsfEleValueMapIsoRhoCut"); + +GsfEleValueMapIsoRhoCut::GsfEleValueMapIsoRhoCut(const edm::ParameterSet& params) : + CutApplicatorWithEventContentBase(params), + slopeTerm_(params,"slopeTerm"), + slopeStart_(params,"slopeStart"), + constTerm_(params,"constTerm"), + rhoEtStart_(params,"rhoEtStart"), + rhoEA_(params,"rhoEA") +{ + auto rho = params.getParameter("rho"); + if(!rho.label().empty()){ + useRho_=true; + contentTags_.emplace("rho",rho); + }else useRho_=false; + + contentTags_.emplace("value",params.getParameter("value")); + +} + +void GsfEleValueMapIsoRhoCut::setConsumes(edm::ConsumesCollector& cc) { + if(useRho_) contentTokens_.emplace("rho",cc.consumes(contentTags_["rho"])); + contentTokens_.emplace("value",cc.consumes >(contentTags_["value"])); +} + +void GsfEleValueMapIsoRhoCut::getEventContent(const edm::EventBase& ev) { + if(useRho_) ev.getByLabel(contentTags_["rho"],rhoHandle_); + ev.getByLabel(contentTags_["value"],valueHandle_); +} + +CutApplicatorBase::result_type +GsfEleValueMapIsoRhoCut:: +operator()(const reco::GsfElectronPtr& cand) const{ + + const float val = (*valueHandle_)[cand]; + + const float et = cand->et(); + const float cutValue = et > slopeStart_(cand) ? slopeTerm_(cand)*(et-slopeStart_(cand)) + constTerm_(cand) : constTerm_(cand); + const float rhoCutValue = getRhoCorr(cand); + + return val < cutValue+rhoCutValue; +} + + +double GsfEleValueMapIsoRhoCut:: +value(const reco::CandidatePtr& cand) const { + reco::GsfElectronPtr ele(cand); + return (*valueHandle_)[cand]; +} + +float GsfEleValueMapIsoRhoCut::getRhoCorr(const reco::GsfElectronPtr& cand)const{ + if(!useRho_) return 0.; + else{ + const double rho = (*rhoHandle_); + return cand->et() >= rhoEtStart_(cand) ? rhoEA_(cand)*rho : 0.; + } +} diff --git a/RecoEgamma/ElectronIdentification/python/Identification/heepElectronID_HEEPV60_80XAOD_cff.py b/RecoEgamma/ElectronIdentification/python/Identification/heepElectronID_HEEPV60_80XAOD_cff.py new file mode 100644 index 0000000000000..ac97d3bf1aed0 --- /dev/null +++ b/RecoEgamma/ElectronIdentification/python/Identification/heepElectronID_HEEPV60_80XAOD_cff.py @@ -0,0 +1,102 @@ +from PhysicsTools.SelectorUtils.centralIDRegistry import central_id_registry + +# Common functions and classes for ID definition are imported here: +from RecoEgamma.ElectronIdentification.Identification.heepElectronID_tools import * + +# +# The HEEP ID cuts V6.0 below are optimized IDs for PHYS14 Scenario PU 20, bx 25ns +# The cut values are taken from the twiki: +# https://twiki.cern.ch/twiki/bin/view/CMS/HEEPElectronIdentificationRun2#Selection_Cuts_HEEP_V5_1 +# (where they may not stay, if a newer version of cuts becomes available for these +# conditions) +# See also the presentation explaining these working points (this will not change): +# [ ... none available for this particular version ... ] + + +# The cut values for the Barrel and Endcap +idName = "heepElectronID-HEEPV60-80XAOD" +WP_HEEP60_EB = HEEP_WorkingPoint_V1( + idName=idName, + dEtaInSeedCut=0.004, + dPhiInCut=0.06, + full5x5SigmaIEtaIEtaCut=9999, + # Two constants for the GsfEleFull5x5E2x5OverE5x5Cut + minE1x5OverE5x5Cut=0.83, + minE2x5OverE5x5Cut=0.94, + # Three constants for the GsfEleHadronicOverEMLinearCut + # cut = constTerm if value < slopeStart + # cut = slopeTerm * (value - slopeStart) + constTerm if value >= slopeStart + hOverESlopeTerm=0.05, + hOverESlopeStart=0.00, + hOverEConstTerm=1.00, + # Three constants for the GsfEleTrkPtIsoCut: + # cut = constTerm if value < slopeStart + # cut = slopeTerm * (value - slopeStart) + constTerm if value >= slopeStart + trkIsoSlopeTerm=0.00, + trkIsoSlopeStart=0.00, + trkIsoConstTerm=5.00, + # Three constants for the GsfEleEmHadD1IsoRhoCut: + # cut = constTerm if value < slopeStart + # cut = slopeTerm * (value - slopeStart) + constTerm if value >= slopeStart + # Also for the same cut, the effective area for the rho correction of the isolation + ehIsoSlopeTerm=0.03, + ehIsoSlopeStart=0.00, + ehIsoConstTerm=2.00, + effAreaForEHIso=0.28, + # other cuts + dxyCut=0.02, + maxMissingHitsCut=1, + ecalDrivenCut=1 + ) + +WP_HEEP60_EE = HEEP_WorkingPoint_V1( + idName=idName, + dEtaInSeedCut=0.006, + dPhiInCut=0.06, + full5x5SigmaIEtaIEtaCut=0.03, + # Two constants for the GsfEleFull5x5E2x5OverE5x5Cut + minE1x5OverE5x5Cut=-1.0, + minE2x5OverE5x5Cut=-1.0, + # Three constants for the GsfEleHadronicOverEMLinearCut + # cut = constTerm if value < slopeStart + # cut = slopeTerm * (value - slopeStart) + constTerm if value >= slopeStart + hOverESlopeTerm=0.05, + hOverESlopeStart=0.00, + hOverEConstTerm=5, + # Three constants for the GsfEleTrkPtIsoCut: + # cut = constTerm if value < slopeStart + # cut = slopeTerm * (value - slopeStart) + constTerm if value >= slopeStart + trkIsoSlopeTerm=0.00, + trkIsoSlopeStart=0.00, + trkIsoConstTerm=5.00, + # Three constants for the GsfEleEmHadD1IsoRhoCut: + # cut = constTerm if value < slopeStart + # cut = slopeTerm * (value - slopeStart) + constTerm if value >= slopeStart + # Also for the same cut, the effective area for the rho correction of the isolation + ehIsoSlopeTerm=0.03, + ehIsoSlopeStart=50.0, + ehIsoConstTerm=2.50, + effAreaForEHIso=0.28, + # other cuts + dxyCut=0.05, + maxMissingHitsCut=1, + ecalDrivenCut=1 + + ) + +# +# Finally, set up VID configuration for all cuts +# +heepElectronID_HEEPV60_80XAOD = configureHEEPElectronID_V60_80XAOD (idName, WP_HEEP60_EB, WP_HEEP60_EE ) + +# +# The MD5 sum numbers below reflect the exact set of cut variables +# and values above. If anything changes, one has to +# 1) comment out the lines below about the registry, +# 2) run "calculateMD5 +# 3) update the MD5 sum strings below and uncomment the lines again. +# + +central_id_registry.register(heepElectronID_HEEPV60_80XAOD.idName,"a7e52edab71dfae8b854b7bce70b56a7") +heepElectronID_HEEPV60_80XAOD.isPOGApproved = cms.untracked.bool(False) + diff --git a/RecoEgamma/ElectronIdentification/python/Identification/heepElectronID_tools.py b/RecoEgamma/ElectronIdentification/python/Identification/heepElectronID_tools.py index 8fabf19234354..92fde5426a75c 100644 --- a/RecoEgamma/ElectronIdentification/python/Identification/heepElectronID_tools.py +++ b/RecoEgamma/ElectronIdentification/python/Identification/heepElectronID_tools.py @@ -192,6 +192,18 @@ def psetGsfEleFull5x5SigmaIEtaIEtaCut(wpEB, wpEE): isIgnored = cms.bool(False) ) +def psetGsfEleFull5x5SigmaIEtaIEtaWithSatCut(wpEB, wpEE): + return cms.PSet( + cutName = cms.string('GsfEleFull5x5SigmaIEtaIEtaWithSatCut'), + maxSigmaIEtaIEtaEB = cms.double( wpEB.full5x5SigmaIEtaIEtaCut ), + maxSigmaIEtaIEtaEE = cms.double( wpEE.full5x5SigmaIEtaIEtaCut ), + maxNrSatCrysIn5x5EB =cms.int32( 0 ), + maxNrSatCrysIn5x5EE =cms.int32( 0 ), + nrSatCrysValueMap = cms.InputTag("heepIDVarValueMaps","eleNrSaturateIn5x5"), + needsAdditionalProducts = cms.bool(True), + + isIgnored = cms.bool(False) + ) # Configure XxX shower shape cuts def psetGsfEleFull5x5E2x5OverE5x5Cut(wpEB, wpEE): return cms.PSet( @@ -205,7 +217,22 @@ def psetGsfEleFull5x5E2x5OverE5x5Cut(wpEB, wpEE): needsAdditionalProducts = cms.bool(False), isIgnored = cms.bool(False) ) - +# Configure XxX shower shape cuts +def psetGsfEleFull5x5E2x5OverE5x5WithSatCut(wpEB, wpEE): + return cms.PSet( + cutName = cms.string('GsfEleFull5x5E2x5OverE5x5WithSatCut'), + # E1x5 / E5x5 + minE1x5OverE5x5EB = cms.double( wpEB.minE1x5OverE5x5Cut ), + minE1x5OverE5x5EE = cms.double( wpEE.minE1x5OverE5x5Cut ), + # E2x5 / E5x5 + minE2x5OverE5x5EB = cms.double( wpEB.minE2x5OverE5x5Cut ), + minE2x5OverE5x5EE = cms.double( wpEE.minE2x5OverE5x5Cut ), + maxNrSatCrysIn5x5EB =cms.int32( 0 ), + maxNrSatCrysIn5x5EE =cms.int32( 0 ), + nrSatCrysValueMap = cms.InputTag("heepIDVarValueMaps","eleNrSaturateIn5x5"), + needsAdditionalProducts = cms.bool(True), + isIgnored = cms.bool(False) + ) # Configure the cut of E/H def psetGsfEleHadronicOverEMLinearCut(wpEB, wpEE) : return cms.PSet( @@ -260,6 +287,28 @@ def psetGsfEleTrkPtIsoRhoCut(wpEB, wpEE): needsAdditionalProducts = cms.bool(True), isIgnored = cms.bool(False) ) +def psetGsfEleTrkPtNoJetCoreIsoCut(wpEB, wpEE): + return cms.PSet( + cutName = cms.string('GsfEleValueMapIsoRhoCut'), + # Three constants for the GsfEleTrkPtIsoCut + # cut = constTerm if value < slopeStart + # cut = slopeTerm * (value - slopeStart) + constTerm if value >= slopeStart + slopeTermEB = cms.double( wpEB.trkIsoSlopeTerm ), + slopeTermEE = cms.double( wpEE.trkIsoSlopeTerm ), + slopeStartEB = cms.double( wpEB.trkIsoSlopeStart ), + slopeStartEE = cms.double( wpEE.trkIsoSlopeStart ), + constTermEB = cms.double( wpEB.trkIsoConstTerm ), + constTermEE = cms.double( wpEE.trkIsoConstTerm ), + #no rho so we zero it out, if the input tag is empty, its ignored anyways + rhoEtStartEB = cms.double( 999999.), + rhoEtStartEE = cms.double( 999999.), + rhoEAEB = cms.double( 0. ), + rhoEAEE = cms.double( 0. ), + rho = cms.InputTag(""), + value = cms.InputTag("heepIDVarValueMaps","eleTrkPtIsoNoJetCore"), + needsAdditionalProducts = cms.bool(True), + isIgnored = cms.bool(False) + ) # Configure the cut on the EM + Had_depth_1 isolation with rho correction def psetGsfEleEmHadD1IsoRhoCut(wpEB, wpEE): @@ -364,6 +413,31 @@ def configureHEEPElectronID_V60(wpEB, wpEE): ) return parameterSet +def configureHEEPElectronID_V60_80XAOD(idName, wpEB, wpEE): + """ + This function configures the full cms.PSet for a VID ID and returns it. + The inputs: two objects of the type HEEP_WorkingPoint_V1, one + containing the cuts for the Barrel (EB) and the other one for the Endcap (EE). + """ + parameterSet = cms.PSet( + idName = cms.string(idName), + cutFlow = cms.VPSet( + psetMinPtCut(), #0 + psetGsfEleSCEtaMultiRangeCut(), #1 + psetGsfEleDEtaInSeedCut(wpEB,wpEE), #2 + psetGsfEleDPhiInCut(wpEB,wpEE), #3 + psetGsfEleFull5x5SigmaIEtaIEtaWithSatCut(wpEB,wpEE), #4 + psetGsfEleFull5x5E2x5OverE5x5WithSatCut(wpEB,wpEE), #5 + psetGsfEleHadronicOverEMLinearCut(wpEB,wpEE), #6 + psetGsfEleTrkPtNoJetCoreIsoCut(wpEB,wpEE), #7 + psetGsfEleEmHadD1IsoRhoCut(wpEB,wpEE), #8 + psetGsfEleDxyCut(wpEB,wpEE), #9 + psetGsfEleMissingHitsCut(wpEB,wpEE), #10, + psetGsfEleEcalDrivenCut(wpEB,wpEE) #11 + ) + ) + return parameterSet + def configureHEEPElectronID_V61(wpEB, wpEE): """ This function configures the full cms.PSet for a VID ID and returns it. From 2cd9f45ed8eae3019a764c49a249ba58684ec3b7 Mon Sep 17 00:00:00 2001 From: Sam Harper Date: Mon, 6 Jun 2016 12:12:26 +0100 Subject: [PATCH 08/18] adding 1st version of HEEP V6.0 in 80X --- .../interface/EBEECutValues.h | 20 ++++++++++++------- 1 file changed, 13 insertions(+), 7 deletions(-) diff --git a/RecoEgamma/ElectronIdentification/interface/EBEECutValues.h b/RecoEgamma/ElectronIdentification/interface/EBEECutValues.h index 39b087015b84a..3fd8ef0215b80 100644 --- a/RecoEgamma/ElectronIdentification/interface/EBEECutValues.h +++ b/RecoEgamma/ElectronIdentification/interface/EBEECutValues.h @@ -9,21 +9,27 @@ namespace reco { typedef edm::Ptr GsfElectronPtr; } -class EBEECutValues { +template +class EBEECutValuesT { private: - double barrel_; - double endcap_; + T barrel_; + T endcap_; const double barrelCutOff_=1.479; //this is currrently used to identify if object is barrel or endcap but may change public: - EBEECutValues(const edm::ParameterSet& params,const std::string& name): - barrel_(params.getParameter(name+"EB")), - endcap_(params.getParameter(name+"EE")){} - double operator()(const reco::GsfElectronPtr& cand)const{return isBarrel(cand) ? barrel_ : endcap_;} + EBEECutValuesT(const edm::ParameterSet& params,const std::string& name): + barrel_(params.getParameter(name+"EB")), + endcap_(params.getParameter(name+"EE")){} + T operator()(const reco::GsfElectronPtr& cand)const{return isBarrel(cand) ? barrel_ : endcap_;} private: const bool isBarrel(const reco::GsfElectronPtr& cand)const{return std::abs(cand->superCluster()->position().eta()) EBEECutValues; +typedef EBEECutValuesT EBEECutValuesInt; + + + #endif From 333f7c96c4030519115756c6ae520d4d5a86fcdf Mon Sep 17 00:00:00 2001 From: Sam Harper Date: Mon, 6 Jun 2016 12:50:06 +0100 Subject: [PATCH 09/18] adding in parameter set descriptions --- .../plugins/ElectronHEEPIDValueMapProducer.cc | 22 ++++++++++++++++--- 1 file changed, 19 insertions(+), 3 deletions(-) diff --git a/RecoEgamma/ElectronIdentification/plugins/ElectronHEEPIDValueMapProducer.cc b/RecoEgamma/ElectronIdentification/plugins/ElectronHEEPIDValueMapProducer.cc index a7028863dced3..5ba87782641f5 100644 --- a/RecoEgamma/ElectronIdentification/plugins/ElectronHEEPIDValueMapProducer.cc +++ b/RecoEgamma/ElectronIdentification/plugins/ElectronHEEPIDValueMapProducer.cc @@ -171,10 +171,26 @@ void ElectronHEEPIDValueMapProducer::writeValueMap(edm::Event &iEvent, } void ElectronHEEPIDValueMapProducer::fillDescriptions(edm::ConfigurationDescriptions& descriptions) { - //The following says we do not know what parameters are allowed so do no validation - // Please change this to state exactly what you do use, even if it is no parameters + edm::ParameterSetDescription desc; - desc.setUnknown(); + desc.add("ebRecHits",edm::InputTag("reducedEcalRecHitsEB")); + desc.add("eeRecHits",edm::InputTag("reducedEcalRecHitsEE")); + desc.add("beamSpot",edm::InputTag("offlineBeamSpot")); + desc.add("tracks",edm::InputTag("generalTracks")); + desc.add("eles",edm::InputTag("gedGsfElectrons")); + + edm::ParameterSetDescription trkIsoDesc; + trkIsoDesc.add("extRadius",0.3); + trkIsoDesc.add("intRadiusBarrel",0.015); + trkIsoDesc.add("intRadiusEndcap",0.015); + trkIsoDesc.add("stripBarrel",0.015); + trkIsoDesc.add("stripEndcap",0.015); + trkIsoDesc.add("ptMin",0.7); + trkIsoDesc.add("maxVtxDist",0.2); + trkIsoDesc.add("drb",999999999.0); + + desc.add("trkIsoConfig",trkIsoDesc); + descriptions.addDefault(desc); } From 3de6992e6b5bef9ac65526b1db6c03762b2fb426 Mon Sep 17 00:00:00 2001 From: Sam Harper Date: Wed, 8 Jun 2016 12:36:33 +0100 Subject: [PATCH 10/18] fixing bug where it the saturation was causing it to auto fail not auto pass --- .../plugins/cuts/GsfEleFull5x5E2x5OverE5x5WithSatCut.cc | 2 +- .../plugins/cuts/GsfEleFull5x5SigmaIEtaIEtaWithSatCut.cc | 2 +- 2 files changed, 2 insertions(+), 2 deletions(-) diff --git a/RecoEgamma/ElectronIdentification/plugins/cuts/GsfEleFull5x5E2x5OverE5x5WithSatCut.cc b/RecoEgamma/ElectronIdentification/plugins/cuts/GsfEleFull5x5E2x5OverE5x5WithSatCut.cc index 2ff6d17f4c258..9fe972080b0cc 100644 --- a/RecoEgamma/ElectronIdentification/plugins/cuts/GsfEleFull5x5E2x5OverE5x5WithSatCut.cc +++ b/RecoEgamma/ElectronIdentification/plugins/cuts/GsfEleFull5x5E2x5OverE5x5WithSatCut.cc @@ -52,7 +52,7 @@ CutApplicatorBase::result_type GsfEleFull5x5E2x5OverE5x5WithSatCut:: operator()(const reco::GsfElectronPtr& cand) const{ - if((*nrSatCrysValueMap_)[cand]>maxNrSatCrysIn5x5Cut_(cand)) return false; + if((*nrSatCrysValueMap_)[cand]>maxNrSatCrysIn5x5Cut_(cand)) return true; const double e5x5 = cand->full5x5_e5x5(); const double e2x5OverE5x5 = e5x5!=0 ? cand->full5x5_e2x5Max()/e5x5 : 0; diff --git a/RecoEgamma/ElectronIdentification/plugins/cuts/GsfEleFull5x5SigmaIEtaIEtaWithSatCut.cc b/RecoEgamma/ElectronIdentification/plugins/cuts/GsfEleFull5x5SigmaIEtaIEtaWithSatCut.cc index 78d197966e905..25b9bad25777e 100644 --- a/RecoEgamma/ElectronIdentification/plugins/cuts/GsfEleFull5x5SigmaIEtaIEtaWithSatCut.cc +++ b/RecoEgamma/ElectronIdentification/plugins/cuts/GsfEleFull5x5SigmaIEtaIEtaWithSatCut.cc @@ -51,7 +51,7 @@ CutApplicatorBase::result_type GsfEleFull5x5SigmaIEtaIEtaWithSatCut:: operator()(const reco::GsfElectronPtr& cand) const{ - if((*nrSatCrysValueMap_)[cand]>maxNrSatCrysIn5x5Cut_(cand)) return false; + if((*nrSatCrysValueMap_)[cand]>maxNrSatCrysIn5x5Cut_(cand)) return true; else return cand->full5x5_sigmaIetaIeta() < maxSigmaIEtaIEtaCut_(cand); } From e8023cd4e5f768ea235d7ab0c255a6dbf2c901b3 Mon Sep 17 00:00:00 2001 From: Sam Harper Date: Fri, 25 Nov 2016 16:24:52 +0000 Subject: [PATCH 11/18] VID Implimentation of HEEP V7.0 --- .../python/tools/vid_id_tools.py | 10 + RecoEgamma/EgammaIsolationAlgos/BuildFile.xml | 1 + .../interface/EleTkIsolFromCands.h | 77 ++++++++ .../src/EleTkIsolFromCands.cc | 160 ++++++++++++++++ .../plugins/ElectronHEEPIDValueMapProducer.cc | 172 +++++++++++------- .../heepElectronID_HEEPV70_cff.py | 102 +++++++++++ .../Identification/heepElectronID_tools.py | 74 +++++++- .../python/heepIdVarValueMapProducer_cfi.py | 42 +++++ 8 files changed, 572 insertions(+), 66 deletions(-) create mode 100644 RecoEgamma/EgammaIsolationAlgos/interface/EleTkIsolFromCands.h create mode 100644 RecoEgamma/EgammaIsolationAlgos/src/EleTkIsolFromCands.cc create mode 100644 RecoEgamma/ElectronIdentification/python/Identification/heepElectronID_HEEPV70_cff.py create mode 100644 RecoEgamma/ElectronIdentification/python/heepIdVarValueMapProducer_cfi.py diff --git a/PhysicsTools/SelectorUtils/python/tools/vid_id_tools.py b/PhysicsTools/SelectorUtils/python/tools/vid_id_tools.py index f838fc93640e5..7059a2194f34c 100644 --- a/PhysicsTools/SelectorUtils/python/tools/vid_id_tools.py +++ b/PhysicsTools/SelectorUtils/python/tools/vid_id_tools.py @@ -96,6 +96,16 @@ def setupVIDElectronSelection(process,cutflow,patProducer=None,addUserData=True) idName = cutflow.idName.value() addVIDSelectionToPATProducer(patProducer,'egmGsfElectronIDs',idName,addUserData) + #we know cutflow has a name otherwise an exception would have been thrown in setupVIDSelection + if cutflow.idName.value().find("HEEP")!=-1: + #not the ideal way but currently the easiest + useMiniAOD = process.egmGsfElectronIDs.physicsObjectSrc == cms.InputTag('slimmedElectrons') + + from RecoEgamma.ElectronIdentification.Identification.heepElectronID_tools import addHEEPProducersToSeq + addHEEPProducersToSeq(process=process,seq=process.egmGsfElectronIDSequence, + insertIndex=process.egmGsfElectronIDSequence.index(process.egmGsfElectronIDs), + useMiniAOD=useMiniAOD) + #### # Muons #### diff --git a/RecoEgamma/EgammaIsolationAlgos/BuildFile.xml b/RecoEgamma/EgammaIsolationAlgos/BuildFile.xml index c92977ecf09ea..254ee40b15074 100644 --- a/RecoEgamma/EgammaIsolationAlgos/BuildFile.xml +++ b/RecoEgamma/EgammaIsolationAlgos/BuildFile.xml @@ -8,6 +8,7 @@ + diff --git a/RecoEgamma/EgammaIsolationAlgos/interface/EleTkIsolFromCands.h b/RecoEgamma/EgammaIsolationAlgos/interface/EleTkIsolFromCands.h new file mode 100644 index 0000000000000..edbc2318772b9 --- /dev/null +++ b/RecoEgamma/EgammaIsolationAlgos/interface/EleTkIsolFromCands.h @@ -0,0 +1,77 @@ +#ifndef RECOEGAMMA_EGAMMAISOLATIONALGOS_ELETKISOLFROMCANDS_H +#define RECOEGAMMA_EGAMMAISOLATIONALGOS_ELETKISOLFROMCANDS_H + +#include "DataFormats/TrackReco/interface/TrackBase.h" +#include "DataFormats/PatCandidates/interface/PackedCandidate.h" +#include "DataFormats/EgammaCandidates/interface/GsfElectron.h" +#include "DataFormats/Common/interface/View.h" +#include "FWCore/ParameterSet/interface/ParameterSet.h" +#include "FWCore/ParameterSet/interface/ParameterSetDescription.h" + +class EleTkIsolFromCands { + +private: + struct TrkCuts { + float minPt; + float minDR2; + float maxDR2; + float minDEta; + float maxDZ; + float minHits; + float minPixelHits; + float maxDPtPt; + std::vector allowedQualities; + std::vector algosToReject; + explicit TrkCuts(const edm::ParameterSet& para); + static edm::ParameterSetDescription pSetDescript(); + }; + + TrkCuts barrelCuts_,endcapCuts_; + +public: + explicit EleTkIsolFromCands(const edm::ParameterSet& para); + EleTkIsolFromCands(const EleTkIsolFromCands&)=default; + ~EleTkIsolFromCands()=default; + EleTkIsolFromCands& operator=(const EleTkIsolFromCands&)=default; + + static edm::ParameterSetDescription pSetDescript(); + + std::pair calIsol(const reco::TrackBase& trk,const pat::PackedCandidateCollection& cands,const edm::View& eles); + + std::pair calIsol(const double eleEta,const double elePhi,const double eleVZ, + const pat::PackedCandidateCollection& cands, + const edm::View& eles); + + double calIsolPt(const reco::TrackBase& trk,const pat::PackedCandidateCollection& cands, + const edm::View& eles){ + return calIsol(trk,cands,eles).second; + } + + double calIsolPt(const double eleEta,const double elePhi,const double eleVZ, + const pat::PackedCandidateCollection& cands, + const edm::View& eles){ + return calIsol(eleEta,elePhi,eleVZ,cands,eles).second; + } + + static bool passTrkSel(const reco::Track& trk, + const double trkPt, + const TrkCuts& cuts, + const double eleEta,const double elePhi, + const double eleVZ); + +private: + //no qualities specified, accept all, ORed + static bool passQual(const reco::TrackBase& trk, + const std::vector& quals); + static bool passAlgo(const reco::TrackBase& trk, + const std::vector& algosToRej); + //for PF electron candidates the "trk pt" is not the track pt + //its the trk-calo gsfele combination energy * trk sin(theta) + //so the idea is to match with the gsf electron and get the orginal + //gsftrack's pt + double getTrkPt(const reco::TrackBase& trk, + const edm::View& eles); +}; + + +#endif diff --git a/RecoEgamma/EgammaIsolationAlgos/src/EleTkIsolFromCands.cc b/RecoEgamma/EgammaIsolationAlgos/src/EleTkIsolFromCands.cc new file mode 100644 index 0000000000000..4d56a8a4294aa --- /dev/null +++ b/RecoEgamma/EgammaIsolationAlgos/src/EleTkIsolFromCands.cc @@ -0,0 +1,160 @@ +#include "RecoEgamma/EgammaIsolationAlgos/interface/EleTkIsolFromCands.h" + +#include "DataFormats/Math/interface/deltaR.h" +#include "DataFormats/GsfTrackReco/interface/GsfTrack.h" + +EleTkIsolFromCands::TrkCuts::TrkCuts(const edm::ParameterSet& para) +{ + minPt = para.getParameter("minPt"); + auto sq = [](double val){return val*val;}; + minDR2 = sq(para.getParameter("minDR")); + maxDR2 = sq(para.getParameter("maxDR")); + minDEta = para.getParameter("minDEta"); + maxDZ = para.getParameter("maxDZ"); + minHits = para.getParameter("minHits"); + minPixelHits = para.getParameter("minPixelHits"); + maxDPtPt = para.getParameter("maxDPtPt"); + + auto qualNames = para.getParameter >("allowedQualities"); + auto algoNames = para.getParameter >("algosToReject"); + + for(auto& qualName : qualNames){ + allowedQualities.push_back(reco::TrackBase::qualityByName(qualName)); + } + for(auto& algoName : algoNames){ + algosToReject.push_back(reco::TrackBase::algoByName(algoName)); + } + std::sort(algosToReject.begin(),algosToReject.end()); + +} + +edm::ParameterSetDescription EleTkIsolFromCands::TrkCuts::pSetDescript() +{ + edm::ParameterSetDescription desc; + desc.add("minPt",1.0); + desc.add("maxDR",0.3); + desc.add("minDR",0.000); + desc.add("minDEta",0.005); + desc.add("maxDZ",0.1); + desc.add("maxDPtPt",-1); + desc.add("minHits",8); + desc.add("minPixelHits",1); + desc.add >("allowedQualities"); + desc.add >("algosToReject"); + return desc; +} + +EleTkIsolFromCands::EleTkIsolFromCands(const edm::ParameterSet& para): + barrelCuts_(para.getParameter("barrelCuts")), + endcapCuts_(para.getParameter("endcapCuts")) +{ + + +} + +edm::ParameterSetDescription EleTkIsolFromCands::pSetDescript() +{ + edm::ParameterSetDescription desc; + desc.add("barrelCuts",TrkCuts::pSetDescript()); + desc.add("endcapCuts",TrkCuts::pSetDescript()); + return desc; +} + +std::pair +EleTkIsolFromCands::calIsol(const reco::TrackBase& eleTrk, + const pat::PackedCandidateCollection& cands, + const edm::View& eles) +{ + return calIsol(eleTrk.eta(),eleTrk.phi(),eleTrk.vz(),cands,eles); +} + +std::pair +EleTkIsolFromCands::calIsol(const double eleEta,const double elePhi, + const double eleVZ, + const pat::PackedCandidateCollection& cands, + const edm::View& eles) +{ + + double ptSum=0.; + int nrTrks=0; + + const TrkCuts& cuts = std::abs(eleEta)<1.5 ? barrelCuts_ : endcapCuts_; + + for(auto& cand : cands){ + if(cand.charge()!=0){ + const reco::Track& trk = cand.pseudoTrack(); + double trkPt = std::abs(cand.pdgId())!=11 ? trk.pt() : getTrkPt(trk,eles); + if(passTrkSel(trk,trkPt,cuts,eleEta,elePhi,eleVZ)){ + ptSum+=trkPt; + nrTrks++; + } + } + } + return {nrTrks,ptSum}; +} + + +bool EleTkIsolFromCands::passTrkSel(const reco::Track& trk, + const double trkPt,const TrkCuts& cuts, + const double eleEta,const double elePhi, + const double eleVZ) +{ + const float dR2 = reco::deltaR2(eleEta,elePhi,trk.eta(),trk.phi()); + const float dEta = trk.eta()-eleEta; + const float dZ = eleVZ - trk.vz(); + + return dR2>=cuts.minDR2 && dR2<=cuts.maxDR2 && + std::abs(dEta)>=cuts.minDEta && + std::abs(dZ)= cuts.minHits && + trk.hitPattern().numberOfValidPixelHits() >=cuts.minPixelHits && + (trk.ptError()/trkPt < cuts.maxDPtPt || cuts.maxDPtPt<0) && + passQual(trk,cuts.allowedQualities) && + passAlgo(trk,cuts.algosToReject) && + trkPt > cuts.minPt; +} + + + +bool EleTkIsolFromCands:: +passQual(const reco::TrackBase& trk, + const std::vector& quals) +{ + if(quals.empty()) return true; + else{ + for(auto qual : quals) { + if(trk.quality(qual)) return true; + } + return false; + } + +} + +bool EleTkIsolFromCands:: +passAlgo(const reco::TrackBase& trk, + const std::vector& algosToRej) +{ + return algosToRej.empty() || !std::binary_search(algosToRej.begin(),algosToRej.end(),trk.algo()); +} + +//so the working theory here is that the track we have is the electrons gsf track +//if so, lets get the pt of the gsf track before E/p combinations +//if no match found to a gsf ele with a gsftrack, return the pt of the input track +double EleTkIsolFromCands:: +getTrkPt(const reco::TrackBase& trk, + const edm::View& eles) +{ + auto match=[](const reco::TrackBase& trk,const reco::GsfElectron& ele){ + return std::abs(trk.eta()-ele.gsfTrack()->eta())<0.001 && + std::abs(trk.phi()-ele.gsfTrack()->phi())<0.001;// && + }; + for(auto& ele : eles){ + if(ele.gsfTrack().isNonnull()){ + if(match(trk,ele)){ + return ele.gsfTrack()->pt(); + } + } + } + return trk.pt(); + +} diff --git a/RecoEgamma/ElectronIdentification/plugins/ElectronHEEPIDValueMapProducer.cc b/RecoEgamma/ElectronIdentification/plugins/ElectronHEEPIDValueMapProducer.cc index 5ba87782641f5..e3c654810296f 100644 --- a/RecoEgamma/ElectronIdentification/plugins/ElectronHEEPIDValueMapProducer.cc +++ b/RecoEgamma/ElectronIdentification/plugins/ElectronHEEPIDValueMapProducer.cc @@ -2,22 +2,26 @@ #include "FWCore/Framework/interface/stream/EDProducer.h" #include "FWCore/Framework/interface/Event.h" +#include "FWCore/Framework/interface/ESHandle.h" #include "FWCore/Framework/interface/MakerMacros.h" #include "FWCore/ParameterSet/interface/ParameterSet.h" #include "DataFormats/Common/interface/ValueMap.h" #include "DataFormats/Common/interface/View.h" +#include "DataFormats/Common/interface/Handle.h" #include "DataFormats/EgammaCandidates/interface/GsfElectron.h" #include "DataFormats/EgammaCandidates/interface/GsfElectronFwd.h" #include "DataFormats/EcalDetId/interface/EcalSubdetector.h" +#include "DataFormats/PatCandidates/interface/PackedCandidate.h" +#include "DataFormats/GsfTrackReco/interface/GsfTrack.h" #include "Geometry/CaloTopology/interface/CaloTopology.h" #include "Geometry/CaloEventSetup/interface/CaloTopologyRecord.h" #include "RecoEcal/EgammaCoreTools/interface/EcalClusterTools.h" -#include "RecoEgamma/EgammaIsolationAlgos/interface/ElectronTkIsolation.h" +#include "RecoEgamma/EgammaIsolationAlgos/interface/EleTkIsolFromCands.h" #include #include @@ -35,7 +39,11 @@ class ElectronHEEPIDValueMapProducer : public edm::stream::EDProducer<> { static void fillDescriptions(edm::ConfigurationDescriptions& descriptions); private: - + template + struct DualToken { + edm::EDGetTokenT aod; + edm::EDGetTokenT miniAOD; + }; virtual void produce(edm::Event&, const edm::EventSetup&) override; template @@ -49,64 +57,93 @@ class ElectronHEEPIDValueMapProducer : public edm::stream::EDProducer<> { edm::Handle& eeHits, edm::ESHandle& caloTopo); + float calTrkIso(const reco::GsfElectron& ele, + const edm::View& eles, + const std::vector >& handles); + template void setToken(edm::EDGetTokenT& token,edm::InputTag tag){token=consumes(tag);} template void setToken(edm::EDGetTokenT& token,const edm::ParameterSet& iPara,const std::string& tag){token=consumes(iPara.getParameter(tag));} + template void setToken(std::vector >& tokens,const edm::ParameterSet& iPara,const std::string& tagName){ + auto tags =iPara.getParameter >(tagName); + for(auto& tag : tags) { + edm::EDGetTokenT token; + setToken(token,tag); + tokens.push_back(token); + } + } + template void setToken(DualToken& token,const edm::ParameterSet& iPara,const std::string& tagAOD,const std::string& tagMiniAOD){ + token.aod=consumes(iPara.getParameter(tagAOD)); + token.miniAOD=consumes(iPara.getParameter(tagMiniAOD)); + } + template void setToken(std::vector >& tokens,const edm::ParameterSet& iPara,const std::string& tagAOD,const std::string& tagMiniAOD){ + auto tagsAOD =iPara.getParameter >(tagAOD); + auto tagsMiniAOD =iPara.getParameter >(tagMiniAOD); + size_t maxSize = std::max(tagsAOD.size(),tagsMiniAOD.size()); + tokens.clear(); + tokens.resize(maxSize); + for(size_t tagNr=0;tagNr edm::Handle getHandle(const edm::Event& iEvent,const edm::EDGetTokenT& token){ edm::Handle handle; iEvent.getByToken(token,handle); return handle; } - + template edm::Handle getHandle(const edm::Event& iEvent,const DualToken& token){ + edm::Handle handle; + iEvent.getByToken(token.aod,handle); + if(!handle.isValid()) iEvent.getByToken(token.miniAOD,handle); + return handle; + } - edm::EDGetTokenT ebRecHitToken_; - edm::EDGetTokenT eeRecHitToken_; - edm::EDGetTokenT > eleToken_; + template std::vector > getHandles(const edm::Event& iEvent,const std::vector >& tokens){ + std::vector > handles(tokens.size()); + if(tokens.empty()) return handles; + if(!tokens[0].aod.isUninitialized()) iEvent.getByToken(tokens[0].aod,handles[0]); + bool isAOD = handles[0].isValid(); + if(!isAOD && !tokens[0].miniAOD.isUninitialized() ) iEvent.getByToken(tokens[0].miniAOD,handles[0]); + + for(size_t tokenNr=1;tokenNr ebRecHitToken_; + DualToken eeRecHitToken_; + DualToken > eleToken_; + std::vector >candTokens_; edm::EDGetTokenT beamSpotToken_; - edm::EDGetTokenT trkToken_; - - struct TrkIsoParam { - double extRadius; - double intRadiusBarrel; - double intRadiusEndcap; - double stripBarrel; - double stripEndcap; - double ptMin; - double maxVtxDist; - double drb; - TrkIsoParam(const edm::ParameterSet& iPara); - }; - const TrkIsoParam trkIsoParam_; - static const std::string eleTrkPtIsoNoJetCoreLabel_; + + EleTkIsolFromCands trkIsoCalc_; + + static const std::string eleTrkPtIsoLabel_; static const std::string eleNrSaturateIn5x5Label_; }; -const std::string ElectronHEEPIDValueMapProducer::eleTrkPtIsoNoJetCoreLabel_="eleTrkPtIsoNoJetCore"; +const std::string ElectronHEEPIDValueMapProducer::eleTrkPtIsoLabel_="eleTrkPtIso"; const std::string ElectronHEEPIDValueMapProducer::eleNrSaturateIn5x5Label_="eleNrSaturateIn5x5"; -ElectronHEEPIDValueMapProducer::TrkIsoParam::TrkIsoParam(const edm::ParameterSet& iParam): - extRadius(iParam.getParameter("extRadius")), - intRadiusBarrel(iParam.getParameter("intRadiusBarrel")), - intRadiusEndcap(iParam.getParameter("intRadiusEndcap")), - stripBarrel(iParam.getParameter("stripBarrel")), - stripEndcap(iParam.getParameter("stripEndcap")), - ptMin(iParam.getParameter("ptMin")), - maxVtxDist(iParam.getParameter("maxVtxDist")), - drb(iParam.getParameter("drb")) -{ -} - ElectronHEEPIDValueMapProducer::ElectronHEEPIDValueMapProducer(const edm::ParameterSet& iConfig): - trkIsoParam_(iConfig.getParameter("trkIsoConfig")) + trkIsoCalc_(iConfig.getParameter("trkIsoConfig")) { - setToken(ebRecHitToken_,iConfig,"ebRecHits"); - setToken(eeRecHitToken_,iConfig,"eeRecHits"); - setToken(eleToken_,iConfig,"eles"); - setToken(trkToken_,iConfig,"tracks"); + setToken(ebRecHitToken_,iConfig,"ebRecHitsAOD","ebRecHitsMiniAOD"); + setToken(eeRecHitToken_,iConfig,"eeRecHitsAOD","eeRecHitsMiniAOD"); + setToken(eleToken_,iConfig,"elesAOD","elesMiniAOD"); + setToken(candTokens_,iConfig,"candsAOD","candsMiniAOD"); setToken(beamSpotToken_,iConfig,"beamSpot"); - produces >(eleTrkPtIsoNoJetCoreLabel_); + produces >(eleTrkPtIsoLabel_); produces >(eleNrSaturateIn5x5Label_); } @@ -118,32 +155,26 @@ ElectronHEEPIDValueMapProducer::~ElectronHEEPIDValueMapProducer() void ElectronHEEPIDValueMapProducer::produce(edm::Event& iEvent, const edm::EventSetup& iSetup) { auto eleHandle = getHandle(iEvent,eleToken_); - auto trkHandle = getHandle(iEvent,trkToken_); auto ebRecHitHandle = getHandle(iEvent,ebRecHitToken_); auto eeRecHitHandle = getHandle(iEvent,eeRecHitToken_); auto beamSpotHandle = getHandle(iEvent,beamSpotToken_); + auto candHandles = getHandles(iEvent,candTokens_); + edm::ESHandle caloTopoHandle; iSetup.get().get(caloTopoHandle); - ElectronTkIsolation isolCorr(trkIsoParam_.extRadius,trkIsoParam_.intRadiusBarrel,trkIsoParam_.intRadiusEndcap, - trkIsoParam_.stripBarrel,trkIsoParam_.stripEndcap,trkIsoParam_.ptMin, - trkIsoParam_.maxVtxDist,trkIsoParam_.drb, - trkHandle.product(),beamSpotHandle->position()); - isolCorr.setAlgosToReject({reco::TrackBase::jetCoreRegionalStep}); - - std::vector eleTrkPtIsoNoJetCore; + std::vector eleTrkPtIso; std::vector eleNrSaturateIn5x5; for(size_t eleNr=0;eleNrsize();eleNr++){ auto elePtr = eleHandle->ptrAt(eleNr); - eleTrkPtIsoNoJetCore.push_back(isolCorr.getPtTracks(&*elePtr)); + eleTrkPtIso.push_back(calTrkIso(*elePtr,*eleHandle,candHandles)); eleNrSaturateIn5x5.push_back(nrSaturatedCrysIn5x5(*elePtr,ebRecHitHandle,eeRecHitHandle,caloTopoHandle)); } - writeValueMap(iEvent,eleHandle,eleTrkPtIsoNoJetCore,eleTrkPtIsoNoJetCoreLabel_); + writeValueMap(iEvent,eleHandle,eleTrkPtIso,eleTrkPtIsoLabel_); writeValueMap(iEvent,eleHandle,eleNrSaturateIn5x5,eleNrSaturateIn5x5Label_); - } int ElectronHEEPIDValueMapProducer::nrSaturatedCrysIn5x5(const reco::GsfElectron& ele, @@ -157,6 +188,22 @@ int ElectronHEEPIDValueMapProducer::nrSaturatedCrysIn5x5(const reco::GsfElectron } +float ElectronHEEPIDValueMapProducer::calTrkIso(const reco::GsfElectron& ele, + const edm::View& eles, + const std::vector >& handles) +{ + if(ele.gsfTrack().isNull()) return std::numeric_limits::max(); + else{ + float trkIso=0.; + for(auto& handle: handles){ + if(handle.isValid()){ + trkIso+= trkIsoCalc_.calIsolPt(*ele.gsfTrack(),*handle,eles); + } + } + return trkIso; + } +} + template void ElectronHEEPIDValueMapProducer::writeValueMap(edm::Event &iEvent, const edm::Handle > & handle, @@ -173,23 +220,18 @@ void ElectronHEEPIDValueMapProducer::writeValueMap(edm::Event &iEvent, void ElectronHEEPIDValueMapProducer::fillDescriptions(edm::ConfigurationDescriptions& descriptions) { edm::ParameterSetDescription desc; - desc.add("ebRecHits",edm::InputTag("reducedEcalRecHitsEB")); - desc.add("eeRecHits",edm::InputTag("reducedEcalRecHitsEE")); desc.add("beamSpot",edm::InputTag("offlineBeamSpot")); - desc.add("tracks",edm::InputTag("generalTracks")); - desc.add("eles",edm::InputTag("gedGsfElectrons")); + desc.add("ebRecHitsAOD",edm::InputTag("reducedEcalRecHitsEB")); + desc.add("eeRecHitsAOD",edm::InputTag("reducedEcalRecHitsEE")); + desc.add >("candsAOD",{edm::InputTag("packedCandidates")}); + desc.add("elesAOD",edm::InputTag("gedGsfElectrons")); - edm::ParameterSetDescription trkIsoDesc; - trkIsoDesc.add("extRadius",0.3); - trkIsoDesc.add("intRadiusBarrel",0.015); - trkIsoDesc.add("intRadiusEndcap",0.015); - trkIsoDesc.add("stripBarrel",0.015); - trkIsoDesc.add("stripEndcap",0.015); - trkIsoDesc.add("ptMin",0.7); - trkIsoDesc.add("maxVtxDist",0.2); - trkIsoDesc.add("drb",999999999.0); + desc.add("ebRecHitsMiniAOD",edm::InputTag("reducedEcalRecHitsEB")); + desc.add("eeRecHitsMiniAOD",edm::InputTag("reducedEcalRecHitsEE")); + desc.add >("candsMiniAOD",{edm::InputTag("packedCandidates")}); + desc.add("elesMiniAOD",edm::InputTag("gedGsfElectrons")); - desc.add("trkIsoConfig",trkIsoDesc); + desc.add("trkIsoConfig",EleTkIsolFromCands::pSetDescript()); descriptions.addDefault(desc); } diff --git a/RecoEgamma/ElectronIdentification/python/Identification/heepElectronID_HEEPV70_cff.py b/RecoEgamma/ElectronIdentification/python/Identification/heepElectronID_HEEPV70_cff.py new file mode 100644 index 0000000000000..c19afb7b05501 --- /dev/null +++ b/RecoEgamma/ElectronIdentification/python/Identification/heepElectronID_HEEPV70_cff.py @@ -0,0 +1,102 @@ +from PhysicsTools.SelectorUtils.centralIDRegistry import central_id_registry + +# Common functions and classes for ID definition are imported here: +from RecoEgamma.ElectronIdentification.Identification.heepElectronID_tools import * + +# +# The HEEP ID cuts V6.0 below are optimized IDs for PHYS14 Scenario PU 20, bx 25ns +# The cut values are taken from the twiki: +# https://twiki.cern.ch/twiki/bin/view/CMS/HEEPElectronIdentificationRun2#Selection_Cuts_HEEP_V5_1 +# (where they may not stay, if a newer version of cuts becomes available for these +# conditions) +# See also the presentation explaining these working points (this will not change): +# [ ... none available for this particular version ... ] + + +# The cut values for the Barrel and Endcap +idName = "heepElectronID-HEEPV70" +WP_HEEP70_EB = HEEP_WorkingPoint_V1( + idName=idName, + dEtaInSeedCut=0.004, + dPhiInCut=0.06, + full5x5SigmaIEtaIEtaCut=9999, + # Two constants for the GsfEleFull5x5E2x5OverE5x5Cut + minE1x5OverE5x5Cut=0.83, + minE2x5OverE5x5Cut=0.94, + # Three constants for the GsfEleHadronicOverEMLinearCut + # cut = constTerm if value < slopeStart + # cut = slopeTerm * (value - slopeStart) + constTerm if value >= slopeStart + hOverESlopeTerm=0.05, + hOverESlopeStart=0.00, + hOverEConstTerm=1.00, + # Three constants for the GsfEleTrkPtIsoCut: + # cut = constTerm if value < slopeStart + # cut = slopeTerm * (value - slopeStart) + constTerm if value >= slopeStart + trkIsoSlopeTerm=0.00, + trkIsoSlopeStart=0.00, + trkIsoConstTerm=5.00, + # Three constants for the GsfEleEmHadD1IsoRhoCut: + # cut = constTerm if value < slopeStart + # cut = slopeTerm * (value - slopeStart) + constTerm if value >= slopeStart + # Also for the same cut, the effective area for the rho correction of the isolation + ehIsoSlopeTerm=0.03, + ehIsoSlopeStart=0.00, + ehIsoConstTerm=2.00, + effAreaForEHIso=0.28, + # other cuts + dxyCut=0.02, + maxMissingHitsCut=1, + ecalDrivenCut=1 + ) + +WP_HEEP70_EE = HEEP_WorkingPoint_V1( + idName=idName, + dEtaInSeedCut=0.006, + dPhiInCut=0.06, + full5x5SigmaIEtaIEtaCut=0.03, + # Two constants for the GsfEleFull5x5E2x5OverE5x5Cut + minE1x5OverE5x5Cut=-1.0, + minE2x5OverE5x5Cut=-1.0, + # Three constants for the GsfEleHadronicOverEMLinearCut + # cut = constTerm if value < slopeStart + # cut = slopeTerm * (value - slopeStart) + constTerm if value >= slopeStart + hOverESlopeTerm=0.05, + hOverESlopeStart=0.00, + hOverEConstTerm=5, + # Three constants for the GsfEleTrkPtIsoCut: + # cut = constTerm if value < slopeStart + # cut = slopeTerm * (value - slopeStart) + constTerm if value >= slopeStart + trkIsoSlopeTerm=0.00, + trkIsoSlopeStart=0.00, + trkIsoConstTerm=5.00, + # Three constants for the GsfEleEmHadD1IsoRhoCut: + # cut = constTerm if value < slopeStart + # cut = slopeTerm * (value - slopeStart) + constTerm if value >= slopeStart + # Also for the same cut, the effective area for the rho correction of the isolation + ehIsoSlopeTerm=0.03, + ehIsoSlopeStart=50.0, + ehIsoConstTerm=2.50, + effAreaForEHIso=0.28, + # other cuts + dxyCut=0.05, + maxMissingHitsCut=1, + ecalDrivenCut=1 + + ) + +# +# Finally, set up VID configuration for all cuts +# +heepElectronID_HEEPV70 = configureHEEPElectronID_V70 (idName, WP_HEEP70_EB, WP_HEEP70_EE ) + +# +# The MD5 sum numbers below reflect the exact set of cut variables +# and values above. If anything changes, one has to +# 1) comment out the lines below about the registry, +# 2) run "calculateMD5 +# 3) update the MD5 sum strings below and uncomment the lines again. +# + +central_id_registry.register(heepElectronID_HEEPV70.idName,"45d49ea5f46f3f13f579d208e5e3412b") +heepElectronID_HEEPV70.isPOGApproved = cms.untracked.bool(False) + diff --git a/RecoEgamma/ElectronIdentification/python/Identification/heepElectronID_tools.py b/RecoEgamma/ElectronIdentification/python/Identification/heepElectronID_tools.py index 92fde5426a75c..a753d74144c8a 100644 --- a/RecoEgamma/ElectronIdentification/python/Identification/heepElectronID_tools.py +++ b/RecoEgamma/ElectronIdentification/python/Identification/heepElectronID_tools.py @@ -309,7 +309,28 @@ def psetGsfEleTrkPtNoJetCoreIsoCut(wpEB, wpEE): needsAdditionalProducts = cms.bool(True), isIgnored = cms.bool(False) ) - +def psetGsfEleTrkPtFall16IsoCut(wpEB, wpEE): + return cms.PSet( + cutName = cms.string('GsfEleValueMapIsoRhoCut'), + # Three constants for the GsfEleTrkPtIsoCut + # cut = constTerm if value < slopeStart + # cut = slopeTerm * (value - slopeStart) + constTerm if value >= slopeStart + slopeTermEB = cms.double( wpEB.trkIsoSlopeTerm ), + slopeTermEE = cms.double( wpEE.trkIsoSlopeTerm ), + slopeStartEB = cms.double( wpEB.trkIsoSlopeStart ), + slopeStartEE = cms.double( wpEE.trkIsoSlopeStart ), + constTermEB = cms.double( wpEB.trkIsoConstTerm ), + constTermEE = cms.double( wpEE.trkIsoConstTerm ), + #no rho so we zero it out, if the input tag is empty, its ignored anyways + rhoEtStartEB = cms.double( 999999.), + rhoEtStartEE = cms.double( 999999.), + rhoEAEB = cms.double( 0. ), + rhoEAEE = cms.double( 0. ), + rho = cms.InputTag(""), + value = cms.InputTag("heepIDVarValueMaps","eleTrkPtIso"), + needsAdditionalProducts = cms.bool(True), + isIgnored = cms.bool(False) + ) # Configure the cut on the EM + Had_depth_1 isolation with rho correction def psetGsfEleEmHadD1IsoRhoCut(wpEB, wpEE): return cms.PSet( @@ -413,6 +434,33 @@ def configureHEEPElectronID_V60(wpEB, wpEE): ) return parameterSet + +def configureHEEPElectronID_V70(idName, wpEB, wpEE): + """ + This function configures the full cms.PSet for a VID ID and returns it. + The inputs: two objects of the type HEEP_WorkingPoint_V1, one + containing the cuts for the Barrel (EB) and the other one for the Endcap (EE). + """ + parameterSet = cms.PSet( + idName = cms.string(idName), + cutFlow = cms.VPSet( + psetMinPtCut(), #0 + psetGsfEleSCEtaMultiRangeCut(), #1 + psetGsfEleDEtaInSeedCut(wpEB,wpEE), #2 + psetGsfEleDPhiInCut(wpEB,wpEE), #3 + psetGsfEleFull5x5SigmaIEtaIEtaWithSatCut(wpEB,wpEE), #4 + psetGsfEleFull5x5E2x5OverE5x5WithSatCut(wpEB,wpEE), #5 + psetGsfEleHadronicOverEMLinearCut(wpEB,wpEE), #6 + psetGsfEleTrkPtFall16IsoCut(wpEB,wpEE), #7 + psetGsfEleEmHadD1IsoRhoCut(wpEB,wpEE), #8 + psetGsfEleDxyCut(wpEB,wpEE), #9 + psetGsfEleMissingHitsCut(wpEB,wpEE), #10, + psetGsfEleEcalDrivenCut(wpEB,wpEE) #11 + ) + ) + return parameterSet + + def configureHEEPElectronID_V60_80XAOD(idName, wpEB, wpEE): """ This function configures the full cms.PSet for a VID ID and returns it. @@ -462,3 +510,27 @@ def configureHEEPElectronID_V61(wpEB, wpEE): ) ) return parameterSet + +def addHEEPProducersToSeq(process,seq,insertIndex,useMiniAOD): + process.load("RecoEgamma.ElectronIdentification.heepIdVarValueMapProducer_cfi") + + seq.insert(insertIndex,process.heepIDVarValueMaps) + + if useMiniAOD==False: + process.load("TrackingTools.TransientTrack.TransientTrackBuilder_cfi") + process.load("PhysicsTools.PatAlgos.slimming.primaryVertexAssociation_cfi") + process.load("PhysicsTools.PatAlgos.slimming.offlineSlimmedPrimaryVertices_cfi") + process.load("PhysicsTools.PatAlgos.slimming.packedPFCandidates_cfi") + from PhysicsTools.PatAlgos.slimming.packedPFCandidates_cfi import packedPFCandidates + process.packedCandsForTkIso = packedPFCandidates.clone() + process.packedCandsForTkIso.PuppiSrc=cms.InputTag("") + process.packedCandsForTkIso.PuppiNoLepSrc=cms.InputTag("") + + process.load("PhysicsTools.PatAlgos.slimming.lostTracks_cfi") + from PhysicsTools.PatAlgos.slimming.lostTracks_cfi import lostTracks + process.lostTracksForTkIso = lostTracks.clone() + process.lostTracksForTkIso.packedPFCandidates =cms.InputTag("packedCandsForTkIso") + seq.insert(insertIndex,process.primaryVertexAssociation) + seq.insert(insertIndex+1,process.offlineSlimmedPrimaryVertices) + seq.insert(insertIndex+2,process.packedCandsForTkIso) + seq.insert(insertIndex+3,process.lostTracksForTkIso) diff --git a/RecoEgamma/ElectronIdentification/python/heepIdVarValueMapProducer_cfi.py b/RecoEgamma/ElectronIdentification/python/heepIdVarValueMapProducer_cfi.py new file mode 100644 index 0000000000000..7d2928875a4ab --- /dev/null +++ b/RecoEgamma/ElectronIdentification/python/heepIdVarValueMapProducer_cfi.py @@ -0,0 +1,42 @@ +import FWCore.ParameterSet.Config as cms + +heepIDVarValueMaps = cms.EDProducer("ElectronHEEPIDValueMapProducer", + beamSpot=cms.InputTag("offlineBeamSpot"), + ebRecHitsAOD=cms.InputTag("reducedEcalRecHitsEB"), + eeRecHitsAOD=cms.InputTag("reducedEcalRecHitsEB"), + candsAOD=cms.VInputTag("packedCandsForTkIso", + "lostTracksForTkIso"), + elesAOD=cms.InputTag("gedGsfElectrons"), + ebRecHitsMiniAOD=cms.InputTag("reducedEgamma","reducedEBRecHits"), + eeRecHitsMiniAOD=cms.InputTag("reducedEgamma","reducedEERecHits"), + candsMiniAOD=cms.VInputTag("packedPFCandidates", + "lostTracks"), + elesMiniAOD=cms.InputTag("slimmedElectrons"), + + trkIsoConfig= cms.PSet( + barrelCuts=cms.PSet( + minPt=cms.double(1.0), + maxDR=cms.double(0.3), + minDR=cms.double(0.0), + minDEta=cms.double(0.005), + maxDZ=cms.double(0.1), + maxDPtPt=cms.double(0.1), + minHits=cms.int32(8), + minPixelHits=cms.int32(1), + allowedQualities=cms.vstring(), + algosToReject=cms.vstring() + ), + endcapCuts=cms.PSet( + minPt=cms.double(1.0), + maxDR=cms.double(0.3), + minDR=cms.double(0.0), + minDEta=cms.double(0.005), + maxDZ=cms.double(0.5), + maxDPtPt=cms.double(0.1), + minHits=cms.int32(8), + minPixelHits=cms.int32(1), + allowedQualities=cms.vstring(), + algosToReject=cms.vstring() + ) + ) + ) From 933ff5b70a421cf581fb88f4a17a00cb3cdb9321 Mon Sep 17 00:00:00 2001 From: Sam Harper Date: Wed, 30 Nov 2016 18:34:35 +0000 Subject: [PATCH 12/18] resetting the ElectronTkIsolation class to its default 80X state --- .../interface/ElectronTkIsolation.h | 11 ++---- .../src/ElectronTkIsolation.cc | 36 +++++++------------ 2 files changed, 16 insertions(+), 31 deletions(-) diff --git a/RecoEgamma/EgammaIsolationAlgos/interface/ElectronTkIsolation.h b/RecoEgamma/EgammaIsolationAlgos/interface/ElectronTkIsolation.h index 350aafd678c10..3add1085d337b 100644 --- a/RecoEgamma/EgammaIsolationAlgos/interface/ElectronTkIsolation.h +++ b/RecoEgamma/EgammaIsolationAlgos/interface/ElectronTkIsolation.h @@ -45,7 +45,7 @@ class ElectronTkIsolation { drb_(drb), trackCollection_(trackCollection), beamPoint_(beamPoint) { - setDefaultAlgosToReject(); + setDzOption("vz"); } @@ -70,7 +70,7 @@ class ElectronTkIsolation { drb_(drb), trackCollection_(trackCollection), beamPoint_(beamPoint) { - setDefaultAlgosToReject(); + setDzOption("vz"); } @@ -104,13 +104,9 @@ class ElectronTkIsolation { double getPtTracks (const reco::GsfElectron*) const ; std::pairgetIso(const reco::GsfElectron*) const; std::pairgetIso(const reco::Track*) const ; - - void setAlgosToReject(std::vector algos); //intentional by value pass private: - - bool passAlgo(const reco::TrackBase& trk)const; - void setDefaultAlgosToReject(){} + double extRadius_ ; double intRadiusBarrel_ ; double intRadiusEndcap_ ; @@ -119,7 +115,6 @@ class ElectronTkIsolation { double ptLow_ ; double lip_ ; double drb_; - std::vector algosToReject_; //vector is sorted const reco::TrackCollection *trackCollection_ ; reco::TrackBase::Point beamPoint_; diff --git a/RecoEgamma/EgammaIsolationAlgos/src/ElectronTkIsolation.cc b/RecoEgamma/EgammaIsolationAlgos/src/ElectronTkIsolation.cc index f52f55e6b1680..92c1eed3fead0 100644 --- a/RecoEgamma/EgammaIsolationAlgos/src/ElectronTkIsolation.cc +++ b/RecoEgamma/EgammaIsolationAlgos/src/ElectronTkIsolation.cc @@ -45,7 +45,6 @@ ElectronTkIsolation::ElectronTkIsolation (double extRadius, trackCollection_(trackCollection), beamPoint_(beamPoint) { - setDefaultAlgosToReject(); setDzOption(dzOptionString); } @@ -87,13 +86,19 @@ std::pair ElectronTkIsolation::getIso(const reco::Track* tmpTrack) c if (fabs( (*itrTr).dxy(beamPoint_) ) > drb_ ) continue; double dr = ROOT::Math::VectorUtil::DeltaR(itrTr->momentum(),tmpElectronMomentumAtVtx) ; double deta = (*itrTr).eta() - tmpElectronEtaAtVertex; - bool isBarrel = std::abs(tmpElectronEtaAtVertex) < 1.479; - double intRadius = isBarrel ? intRadiusBarrel_ : intRadiusEndcap_; - double strip = isBarrel ? stripBarrel_ : stripEndcap_; - if(dr < extRadius_ && dr>=intRadius && std::abs(deta) >=strip && passAlgo(*itrTr)){ - - ++counter ; - ptSum += this_pt; + if (fabs(tmpElectronEtaAtVertex) < 1.479) { + if ( fabs(dr) < extRadius_ && fabs(dr) >= intRadiusBarrel_ && fabs(deta) >= stripBarrel_) + { + ++counter ; + ptSum += this_pt; + } + } + else { + if ( fabs(dr) < extRadius_ && fabs(dr) >= intRadiusEndcap_ && fabs(deta) >= stripEndcap_) + { + ++counter ; + ptSum += this_pt; + } } }//end loop over tracks @@ -117,18 +122,3 @@ double ElectronTkIsolation::getPtTracks (const reco::GsfElectron* electron) cons return getIso(electron).second ; } - -bool ElectronTkIsolation::passAlgo(const reco::TrackBase& trk)const -{ - int algo = trk.algo(); - bool rejAlgo=std::binary_search(algosToReject_.begin(),algosToReject_.end(),algo); - return rejAlgo==false; -} - - -void ElectronTkIsolation::setAlgosToReject(std::vector algos) -{ - algosToReject_=std::move(algos); - std::sort(algosToReject_.begin(),algosToReject_.end()); -} - From fbf760eedf7286f89efbfc34af53d446aa6ae090 Mon Sep 17 00:00:00 2001 From: Sam Harper Date: Thu, 8 Dec 2016 23:50:32 +0000 Subject: [PATCH 13/18] minor fixes for code review --- .../SelectorUtils/python/tools/vid_id_tools.py | 6 +++++- .../src/EleTkIsolFromCands.cc | 15 ++++++++------- .../plugins/cuts/GsfEleE2x5OverE5x5Cut.cc | 6 +++--- .../plugins/cuts/GsfEleFull5x5E2x5OverE5x5Cut.cc | 4 +++- .../cuts/GsfEleFull5x5E2x5OverE5x5WithSatCut.cc | 4 +++- 5 files changed, 22 insertions(+), 13 deletions(-) diff --git a/PhysicsTools/SelectorUtils/python/tools/vid_id_tools.py b/PhysicsTools/SelectorUtils/python/tools/vid_id_tools.py index 7059a2194f34c..f7067152588ab 100644 --- a/PhysicsTools/SelectorUtils/python/tools/vid_id_tools.py +++ b/PhysicsTools/SelectorUtils/python/tools/vid_id_tools.py @@ -97,7 +97,11 @@ def setupVIDElectronSelection(process,cutflow,patProducer=None,addUserData=True) addVIDSelectionToPATProducer(patProducer,'egmGsfElectronIDs',idName,addUserData) #we know cutflow has a name otherwise an exception would have been thrown in setupVIDSelection - if cutflow.idName.value().find("HEEP")!=-1: + #run this for all heep IDs except V60 standard for which it is not needed + #it is needed for V61 and V70 + if (cutflow.idName.value().find("HEEP")!=-1 and + cutflow.idName.value()!="heepElectronID-HEEPV60"): + #not the ideal way but currently the easiest useMiniAOD = process.egmGsfElectronIDs.physicsObjectSrc == cms.InputTag('slimmedElectrons') diff --git a/RecoEgamma/EgammaIsolationAlgos/src/EleTkIsolFromCands.cc b/RecoEgamma/EgammaIsolationAlgos/src/EleTkIsolFromCands.cc index 4d56a8a4294aa..89b73997fcee5 100644 --- a/RecoEgamma/EgammaIsolationAlgos/src/EleTkIsolFromCands.cc +++ b/RecoEgamma/EgammaIsolationAlgos/src/EleTkIsolFromCands.cc @@ -121,13 +121,12 @@ passQual(const reco::TrackBase& trk, const std::vector& quals) { if(quals.empty()) return true; - else{ - for(auto qual : quals) { - if(trk.quality(qual)) return true; - } - return false; + + for(auto qual : quals) { + if(trk.quality(qual)) return true; } - + + return false; } bool EleTkIsolFromCands:: @@ -144,9 +143,11 @@ double EleTkIsolFromCands:: getTrkPt(const reco::TrackBase& trk, const edm::View& eles) { + //note, the trk.eta(),trk.phi() should be identical to the gsf track eta,phi + //although this may not be the case due to roundings after packing auto match=[](const reco::TrackBase& trk,const reco::GsfElectron& ele){ return std::abs(trk.eta()-ele.gsfTrack()->eta())<0.001 && - std::abs(trk.phi()-ele.gsfTrack()->phi())<0.001;// && + reco::deltaPhi(trk.phi(),ele.gsfTrack()->phi())<0.001;// && }; for(auto& ele : eles){ if(ele.gsfTrack().isNonnull()){ diff --git a/RecoEgamma/ElectronIdentification/plugins/cuts/GsfEleE2x5OverE5x5Cut.cc b/RecoEgamma/ElectronIdentification/plugins/cuts/GsfEleE2x5OverE5x5Cut.cc index 0908a86409f5a..1a887fb0064bd 100644 --- a/RecoEgamma/ElectronIdentification/plugins/cuts/GsfEleE2x5OverE5x5Cut.cc +++ b/RecoEgamma/ElectronIdentification/plugins/cuts/GsfEleE2x5OverE5x5Cut.cc @@ -39,7 +39,7 @@ operator()(const reco::GsfElectronPtr& cand) const{ double GsfEleE2x5OverE5x5Cut::value(const reco::CandidatePtr& cand) const { reco::GsfElectronPtr ele(cand); - // here return the ratio since it is a decent encoding of the value of - // both variables - return ele->e2x5Max()/ele->e1x5(); + //btw we broke somebodies nice model of assuming every cut is 1D.... + //what this is returning is fairly meaningless... + return ele->e1x5() ? ele->e2x5Max()/ele->e1x5() : 0.; } diff --git a/RecoEgamma/ElectronIdentification/plugins/cuts/GsfEleFull5x5E2x5OverE5x5Cut.cc b/RecoEgamma/ElectronIdentification/plugins/cuts/GsfEleFull5x5E2x5OverE5x5Cut.cc index 2fc0ab7231059..0e90bf33e936e 100644 --- a/RecoEgamma/ElectronIdentification/plugins/cuts/GsfEleFull5x5E2x5OverE5x5Cut.cc +++ b/RecoEgamma/ElectronIdentification/plugins/cuts/GsfEleFull5x5E2x5OverE5x5Cut.cc @@ -49,5 +49,7 @@ operator()(const reco::GsfElectronPtr& cand) const{ double GsfEleFull5x5E2x5OverE5x5Cut:: value(const reco::CandidatePtr& cand) const { reco::GsfElectronPtr ele(cand); - return ele->full5x5_e2x5Max()/ele->full5x5_e1x5(); + //btw we broke somebodies nice model of assuming every cut is 1D.... + //what this is returning is fairly meaningless... + return ele->full5x5_e1x5() ? ele->full5x5_e2x5Max()/ele->full5x5_e1x5() : 0.; } diff --git a/RecoEgamma/ElectronIdentification/plugins/cuts/GsfEleFull5x5E2x5OverE5x5WithSatCut.cc b/RecoEgamma/ElectronIdentification/plugins/cuts/GsfEleFull5x5E2x5OverE5x5WithSatCut.cc index 9fe972080b0cc..104ec9f9ac55d 100644 --- a/RecoEgamma/ElectronIdentification/plugins/cuts/GsfEleFull5x5E2x5OverE5x5WithSatCut.cc +++ b/RecoEgamma/ElectronIdentification/plugins/cuts/GsfEleFull5x5E2x5OverE5x5WithSatCut.cc @@ -65,5 +65,7 @@ operator()(const reco::GsfElectronPtr& cand) const{ double GsfEleFull5x5E2x5OverE5x5WithSatCut:: value(const reco::CandidatePtr& cand) const { reco::GsfElectronPtr ele(cand); - return ele->full5x5_e2x5Max()/ele->full5x5_e1x5(); + //btw we broke somebodies nice model of assuming every cut is 1D.... + //what this is returning is fairly meaningless... + return ele->full5x5_e1x5() ? ele->full5x5_e2x5Max()/ele->full5x5_e1x5() : 0.; } From abbbb025c9745f699111f94743783414b752aa60 Mon Sep 17 00:00:00 2001 From: Sam Harper Date: Sun, 11 Dec 2016 19:44:45 +0000 Subject: [PATCH 14/18] updating unit tests and switching the pog approved bool to true --- .../python/Identification/heepElectronID_HEEPV70_cff.py | 2 +- RecoEgamma/ElectronIdentification/test/runElectron_VID.py | 3 ++- RecoEgamma/ElectronIdentification/test/runtests.sh | 1 + 3 files changed, 4 insertions(+), 2 deletions(-) diff --git a/RecoEgamma/ElectronIdentification/python/Identification/heepElectronID_HEEPV70_cff.py b/RecoEgamma/ElectronIdentification/python/Identification/heepElectronID_HEEPV70_cff.py index c19afb7b05501..097219cf273e7 100644 --- a/RecoEgamma/ElectronIdentification/python/Identification/heepElectronID_HEEPV70_cff.py +++ b/RecoEgamma/ElectronIdentification/python/Identification/heepElectronID_HEEPV70_cff.py @@ -98,5 +98,5 @@ # central_id_registry.register(heepElectronID_HEEPV70.idName,"45d49ea5f46f3f13f579d208e5e3412b") -heepElectronID_HEEPV70.isPOGApproved = cms.untracked.bool(False) +heepElectronID_HEEPV70.isPOGApproved = cms.untracked.bool(True) diff --git a/RecoEgamma/ElectronIdentification/test/runElectron_VID.py b/RecoEgamma/ElectronIdentification/test/runElectron_VID.py index 3d12ffadd1289..41f9dfb679a22 100644 --- a/RecoEgamma/ElectronIdentification/test/runElectron_VID.py +++ b/RecoEgamma/ElectronIdentification/test/runElectron_VID.py @@ -1,3 +1,4 @@ + import FWCore.ParameterSet.Config as cms import sys @@ -7,7 +8,7 @@ process.MessageLogger.cerr.FwkReport.reportEvery = 100 process.load("Configuration.StandardSequences.GeometryDB_cff") - +process.load("Configuration.StandardSequences.MagneticField_cff") process.load("Configuration.StandardSequences.FrontierConditions_GlobalTag_cff") # NOTE: the pick the right global tag! # for PHYS14 scenario PU4bx50 : global tag is ??? diff --git a/RecoEgamma/ElectronIdentification/test/runtests.sh b/RecoEgamma/ElectronIdentification/test/runtests.sh index b338f0294970f..2648092eb8583 100755 --- a/RecoEgamma/ElectronIdentification/test/runtests.sh +++ b/RecoEgamma/ElectronIdentification/test/runtests.sh @@ -5,6 +5,7 @@ ids_to_test=( "RecoEgamma.ElectronIdentification.Identification.cutBasedElectronID_Spring15_25ns_V1_cff" "RecoEgamma.ElectronIdentification.Identification.cutBasedElectronID_Spring15_50ns_V2_cff" "RecoEgamma.ElectronIdentification.Identification.heepElectronID_HEEPV60_cff" +"RecoEgamma.ElectronIdentification.Identification.heepElectronID_HEEPV70_cff" "RecoEgamma.ElectronIdentification.Identification.mvaElectronID_Spring15_25ns_nonTrig_V1_cff" "RecoEgamma.ElectronIdentification.Identification.mvaElectronID_Spring15_25ns_Trig_V1_cff" "RecoEgamma.ElectronIdentification.Identification.mvaElectronID_Spring15_50ns_Trig_V1_cff" From 2cc5e898904cd1685af8ecd5250e9633ab16a61e Mon Sep 17 00:00:00 2001 From: Sam Harper Date: Mon, 12 Dec 2016 23:23:49 +0000 Subject: [PATCH 15/18] updating to allow HEEP V7.0 to be packed into miniAOD if it is so wished --- PhysicsTools/PatAlgos/python/slimming/miniAOD_tools.py | 7 +++++++ 1 file changed, 7 insertions(+) diff --git a/PhysicsTools/PatAlgos/python/slimming/miniAOD_tools.py b/PhysicsTools/PatAlgos/python/slimming/miniAOD_tools.py index 5fce939a436dc..8ef35407e7174 100644 --- a/PhysicsTools/PatAlgos/python/slimming/miniAOD_tools.py +++ b/PhysicsTools/PatAlgos/python/slimming/miniAOD_tools.py @@ -186,6 +186,13 @@ def miniAOD_customizeCommon(process): cms.InputTag('reducedEgamma','reducedGedGsfElectrons') for idmod in electron_ids: setupAllVIDIdsInModule(process,idmod,setupVIDElectronSelection,None,False) + + #heepIDVarValueMaps only exists if HEEP V6.1 or HEEP 7.0 ID has already been loaded + if hasattr(process,'heepIDVarValueMaps'): + process.heepIDVarValueMaps.elesMiniAOD = cms.InputTag('reducedEgamma','reducedGedGsfElectrons') + #force HEEP to use miniAOD (otherwise it'll detect the AOD) + process.heepIDVarValueMaps.dataFormat = cms.int32(2) + #VID Photon IDs photon_ids = ['RecoEgamma.PhotonIdentification.Identification.cutBasedPhotonID_Spring15_25ns_V1_cff', From 9194d86981afbe6ac4f5d2e8348ce1d645bf00e7 Mon Sep 17 00:00:00 2001 From: Sam Harper Date: Mon, 12 Dec 2016 18:50:03 +0000 Subject: [PATCH 16/18] updating unit tests to a 80X sample (necessary for HEEP ID V70 which needs 80X products --- .../test/runElectron_VID.py | 31 ++++++++++--------- 1 file changed, 16 insertions(+), 15 deletions(-) diff --git a/RecoEgamma/ElectronIdentification/test/runElectron_VID.py b/RecoEgamma/ElectronIdentification/test/runElectron_VID.py index 41f9dfb679a22..d1ed67b37a294 100644 --- a/RecoEgamma/ElectronIdentification/test/runElectron_VID.py +++ b/RecoEgamma/ElectronIdentification/test/runElectron_VID.py @@ -22,23 +22,24 @@ # Define input data to read # process.maxEvents = cms.untracked.PSet( input = cms.untracked.int32(1000) ) - +#AOD test files dataset :/RelValZEE_13/CMSSW_8_0_21-PU25ns_80X_mcRun2_asymptotic_2016_TrancheIV_v6_Tr4GT_v6-v1/GEN-SIM-RECO inputFilesAOD = cms.untracked.vstring( - # AOD test files from /store/relval/CMSSW_7_6_0_pre4/RelValZEE_13/GEN-SIM-RECO/PU25ns_76X_mcRun2_asymptotic_v1-v1 - '/store/relval/CMSSW_7_6_0_pre4/RelValZEE_13/GEN-SIM-RECO/PU25ns_76X_mcRun2_asymptotic_v1-v1/00000/3ADB5D32-DD4F-E511-AC01-002618943811.root', - '/store/relval/CMSSW_7_6_0_pre4/RelValZEE_13/GEN-SIM-RECO/PU25ns_76X_mcRun2_asymptotic_v1-v1/00000/54B6CF34-DD4F-E511-9629-002590596490.root', - '/store/relval/CMSSW_7_6_0_pre4/RelValZEE_13/GEN-SIM-RECO/PU25ns_76X_mcRun2_asymptotic_v1-v1/00000/8043D96A-6C4F-E511-81E7-003048FFD736.root', - '/store/relval/CMSSW_7_6_0_pre4/RelValZEE_13/GEN-SIM-RECO/PU25ns_76X_mcRun2_asymptotic_v1-v1/00000/8E554BD2-6D4F-E511-BFD2-0025905A60DE.root', - '/store/relval/CMSSW_7_6_0_pre4/RelValZEE_13/GEN-SIM-RECO/PU25ns_76X_mcRun2_asymptotic_v1-v1/00000/98EB5C3F-6D4F-E511-910B-0025905A6056.root', - '/store/relval/CMSSW_7_6_0_pre4/RelValZEE_13/GEN-SIM-RECO/PU25ns_76X_mcRun2_asymptotic_v1-v1/00000/9C8CF66A-6C4F-E511-BD02-00259059391E.root', - '/store/relval/CMSSW_7_6_0_pre4/RelValZEE_13/GEN-SIM-RECO/PU25ns_76X_mcRun2_asymptotic_v1-v1/00000/D015FB85-6C4F-E511-88FE-002618943902.root', - '/store/relval/CMSSW_7_6_0_pre4/RelValZEE_13/GEN-SIM-RECO/PU25ns_76X_mcRun2_asymptotic_v1-v1/00000/D873CC62-6C4F-E511-ABBA-0025905B855E.root' - ) + '/store/relval/CMSSW_8_0_21/RelValZEE_13/GEN-SIM-RECO/PU25ns_80X_mcRun2_asymptotic_2016_TrancheIV_v6_Tr4GT_v6-v1/10000/105EBFE7-9198-E611-8604-0025905B85D2.root', +'/store/relval/CMSSW_8_0_21/RelValZEE_13/GEN-SIM-RECO/PU25ns_80X_mcRun2_asymptotic_2016_TrancheIV_v6_Tr4GT_v6-v1/10000/483FD411-9498-E611-9427-0025905B85D2.root', +'/store/relval/CMSSW_8_0_21/RelValZEE_13/GEN-SIM-RECO/PU25ns_80X_mcRun2_asymptotic_2016_TrancheIV_v6_Tr4GT_v6-v1/10000/4ACA8789-A498-E611-A54C-0025905B858E.root', +'/store/relval/CMSSW_8_0_21/RelValZEE_13/GEN-SIM-RECO/PU25ns_80X_mcRun2_asymptotic_2016_TrancheIV_v6_Tr4GT_v6-v1/10000/5A20FA98-9498-E611-9DB0-0CC47A78A33E.root', +'/store/relval/CMSSW_8_0_21/RelValZEE_13/GEN-SIM-RECO/PU25ns_80X_mcRun2_asymptotic_2016_TrancheIV_v6_Tr4GT_v6-v1/10000/6694D992-9098-E611-8461-0CC47A4D768E.root', +'/store/relval/CMSSW_8_0_21/RelValZEE_13/GEN-SIM-RECO/PU25ns_80X_mcRun2_asymptotic_2016_TrancheIV_v6_Tr4GT_v6-v1/10000/7887DD0D-9498-E611-8D98-0CC47A4C8E64.root', +'/store/relval/CMSSW_8_0_21/RelValZEE_13/GEN-SIM-RECO/PU25ns_80X_mcRun2_asymptotic_2016_TrancheIV_v6_Tr4GT_v6-v1/10000/845B6ABD-9098-E611-A2E3-0CC47A78A4BA.root', +'/store/relval/CMSSW_8_0_21/RelValZEE_13/GEN-SIM-RECO/PU25ns_80X_mcRun2_asymptotic_2016_TrancheIV_v6_Tr4GT_v6-v1/10000/A8877C9C-A498-E611-8EFF-0025905A607E.root', -inputFilesMiniAOD = cms.untracked.vstring( - # MiniAOD test files from /store/relval/CMSSW_7_6_0_pre4/RelValZEE_13/MINIAODSIM/PU25ns_76X_mcRun2_asymptotic_v1-v1 - '/store/relval/CMSSW_7_6_0_pre4/RelValZEE_13/MINIAODSIM/PU25ns_76X_mcRun2_asymptotic_v1-v1/00000/BE21962F-DD4F-E511-B681-002354EF3BDF.root', - '/store/relval/CMSSW_7_6_0_pre4/RelValZEE_13/MINIAODSIM/PU25ns_76X_mcRun2_asymptotic_v1-v1/00000/D2B5E032-DD4F-E511-96A4-0025905A610C.root' + + #'file:/opt/ppd/scratch/harper/mcTestFiles/ZToEE_NNPDF30_13TeV-powheg_M_200_400_PUSpring16_80X_mcRun2_asymptotic_2016_v3-v2_FE668C0A-8B09-E611-ACDE-B083FED76DBD.root' + ) +#AOD test files dataset :/RelValZEE_13/CMSSW_8_0_21-PU25ns_80X_mcRun2_asymptotic_2016_TrancheIV_v6_Tr4GT_v6-v1/MINIAODSIM +inputFilesMiniAOD = cms.untracked.vstring( + '/store/relval/CMSSW_8_0_21/RelValZEE_13/MINIAODSIM/PU25ns_80X_mcRun2_asymptotic_2016_TrancheIV_v6_Tr4GT_v6-v1/10000/0E6B945E-A498-E611-BD29-0CC47A78A30E.root', + '/store/relval/CMSSW_8_0_21/RelValZEE_13/MINIAODSIM/PU25ns_80X_mcRun2_asymptotic_2016_TrancheIV_v6_Tr4GT_v6-v1/10000/EE70625F-A498-E611-8CC7-0CC47A4D766C.root', ) # Set up input/output depending on the format From eb319604a5d26221276994cf7e699ffd8c003c02 Mon Sep 17 00:00:00 2001 From: Sam Harper Date: Mon, 12 Dec 2016 22:45:48 +0000 Subject: [PATCH 17/18] replacing the import * with specific modules --- .../python/Identification/heepElectronID_HEEPV60_80XAOD_cff.py | 3 ++- .../python/Identification/heepElectronID_HEEPV70_cff.py | 3 ++- 2 files changed, 4 insertions(+), 2 deletions(-) diff --git a/RecoEgamma/ElectronIdentification/python/Identification/heepElectronID_HEEPV60_80XAOD_cff.py b/RecoEgamma/ElectronIdentification/python/Identification/heepElectronID_HEEPV60_80XAOD_cff.py index ac97d3bf1aed0..32a9b8461c981 100644 --- a/RecoEgamma/ElectronIdentification/python/Identification/heepElectronID_HEEPV60_80XAOD_cff.py +++ b/RecoEgamma/ElectronIdentification/python/Identification/heepElectronID_HEEPV60_80XAOD_cff.py @@ -1,7 +1,8 @@ +import FWCore.ParameterSet.Config as cms from PhysicsTools.SelectorUtils.centralIDRegistry import central_id_registry # Common functions and classes for ID definition are imported here: -from RecoEgamma.ElectronIdentification.Identification.heepElectronID_tools import * +from RecoEgamma.ElectronIdentification.Identification.heepElectronID_tools import HEEP_WorkingPoint_V1,configureHEEPElectronID_V60_80XAOD # # The HEEP ID cuts V6.0 below are optimized IDs for PHYS14 Scenario PU 20, bx 25ns diff --git a/RecoEgamma/ElectronIdentification/python/Identification/heepElectronID_HEEPV70_cff.py b/RecoEgamma/ElectronIdentification/python/Identification/heepElectronID_HEEPV70_cff.py index 097219cf273e7..40875681bcfea 100644 --- a/RecoEgamma/ElectronIdentification/python/Identification/heepElectronID_HEEPV70_cff.py +++ b/RecoEgamma/ElectronIdentification/python/Identification/heepElectronID_HEEPV70_cff.py @@ -1,7 +1,8 @@ +import FWCore.ParameterSet.Config as cms from PhysicsTools.SelectorUtils.centralIDRegistry import central_id_registry # Common functions and classes for ID definition are imported here: -from RecoEgamma.ElectronIdentification.Identification.heepElectronID_tools import * +from RecoEgamma.ElectronIdentification.Identification.heepElectronID_tools import HEEP_WorkingPoint_V1,configureHEEPElectronID_V70 # # The HEEP ID cuts V6.0 below are optimized IDs for PHYS14 Scenario PU 20, bx 25ns From f11ac4a3601d7c0dc616a260643181967bf41a13 Mon Sep 17 00:00:00 2001 From: Sam Harper Date: Mon, 12 Dec 2016 22:47:42 +0000 Subject: [PATCH 18/18] adding option to force AOD or miniAOD --- .../plugins/ElectronHEEPIDValueMapProducer.cc | 67 ++++++++++++------- .../python/heepIdVarValueMapProducer_cfi.py | 1 + 2 files changed, 43 insertions(+), 25 deletions(-) diff --git a/RecoEgamma/ElectronIdentification/plugins/ElectronHEEPIDValueMapProducer.cc b/RecoEgamma/ElectronIdentification/plugins/ElectronHEEPIDValueMapProducer.cc index e3c654810296f..b72c11446c46d 100644 --- a/RecoEgamma/ElectronIdentification/plugins/ElectronHEEPIDValueMapProducer.cc +++ b/RecoEgamma/ElectronIdentification/plugins/ElectronHEEPIDValueMapProducer.cc @@ -30,20 +30,31 @@ class ElectronHEEPIDValueMapProducer : public edm::stream::EDProducer<> { - +private: + //helper classes to handle AOD vs MiniAOD + template + struct DualToken { + edm::EDGetTokenT aod; + edm::EDGetTokenT miniAOD; + }; + class DataFormat { public: - + enum Format{AUTO=0,AOD=1,MINIAOD=2}; + private: + int data_; + public: + DataFormat(int val):data_(val){} + bool tryAOD()const{return data_==AUTO || data_==AOD;} + bool tryMiniAOD()const{return data_==AUTO || data_==MINIAOD;} + int operator()()const{return data_;} + }; +public: explicit ElectronHEEPIDValueMapProducer(const edm::ParameterSet&); ~ElectronHEEPIDValueMapProducer(); static void fillDescriptions(edm::ConfigurationDescriptions& descriptions); - private: - template - struct DualToken { - edm::EDGetTokenT aod; - edm::EDGetTokenT miniAOD; - }; +private: virtual void produce(edm::Event&, const edm::EventSetup&) override; template @@ -71,21 +82,25 @@ class ElectronHEEPIDValueMapProducer : public edm::stream::EDProducer<> { tokens.push_back(token); } } - template void setToken(DualToken& token,const edm::ParameterSet& iPara,const std::string& tagAOD,const std::string& tagMiniAOD){ - token.aod=consumes(iPara.getParameter(tagAOD)); - token.miniAOD=consumes(iPara.getParameter(tagMiniAOD)); + template void setToken(DualToken& token,const edm::ParameterSet& iPara,const std::string& tagAOD,const std::string& tagMiniAOD,DataFormat format){ + if(format.tryAOD()) token.aod=consumes(iPara.getParameter(tagAOD)); + if(format.tryMiniAOD()) token.miniAOD=consumes(iPara.getParameter(tagMiniAOD)); } - template void setToken(std::vector >& tokens,const edm::ParameterSet& iPara,const std::string& tagAOD,const std::string& tagMiniAOD){ + template void setToken(std::vector >& tokens,const edm::ParameterSet& iPara,const std::string& tagAOD,const std::string& tagMiniAOD,DataFormat format){ auto tagsAOD =iPara.getParameter >(tagAOD); auto tagsMiniAOD =iPara.getParameter >(tagMiniAOD); size_t maxSize = std::max(tagsAOD.size(),tagsMiniAOD.size()); tokens.clear(); - tokens.resize(maxSize); - for(size_t tagNr=0;tagNr { } template edm::Handle getHandle(const edm::Event& iEvent,const DualToken& token){ edm::Handle handle; - iEvent.getByToken(token.aod,handle); - if(!handle.isValid()) iEvent.getByToken(token.miniAOD,handle); + if(!token.aod.isUninitialized()) iEvent.getByToken(token.aod,handle); + if(!handle.isValid() && !token.miniAOD.isUninitialized()) iEvent.getByToken(token.miniAOD,handle); return handle; } @@ -124,6 +139,7 @@ class ElectronHEEPIDValueMapProducer : public edm::stream::EDProducer<> { edm::EDGetTokenT beamSpotToken_; EleTkIsolFromCands trkIsoCalc_; + DataFormat dataFormat_; static const std::string eleTrkPtIsoLabel_; static const std::string eleNrSaturateIn5x5Label_; @@ -135,12 +151,13 @@ const std::string ElectronHEEPIDValueMapProducer::eleNrSaturateIn5x5Label_="eleN ElectronHEEPIDValueMapProducer::ElectronHEEPIDValueMapProducer(const edm::ParameterSet& iConfig): - trkIsoCalc_(iConfig.getParameter("trkIsoConfig")) + trkIsoCalc_(iConfig.getParameter("trkIsoConfig")), + dataFormat_(iConfig.getParameter("dataFormat")) { - setToken(ebRecHitToken_,iConfig,"ebRecHitsAOD","ebRecHitsMiniAOD"); - setToken(eeRecHitToken_,iConfig,"eeRecHitsAOD","eeRecHitsMiniAOD"); - setToken(eleToken_,iConfig,"elesAOD","elesMiniAOD"); - setToken(candTokens_,iConfig,"candsAOD","candsMiniAOD"); + setToken(ebRecHitToken_,iConfig,"ebRecHitsAOD","ebRecHitsMiniAOD",dataFormat_); + setToken(eeRecHitToken_,iConfig,"eeRecHitsAOD","eeRecHitsMiniAOD",dataFormat_); + setToken(eleToken_,iConfig,"elesAOD","elesMiniAOD",dataFormat_); + setToken(candTokens_,iConfig,"candsAOD","candsMiniAOD",dataFormat_); setToken(beamSpotToken_,iConfig,"beamSpot"); produces >(eleTrkPtIsoLabel_); @@ -166,7 +183,6 @@ void ElectronHEEPIDValueMapProducer::produce(edm::Event& iEvent, const edm::Even std::vector eleTrkPtIso; std::vector eleNrSaturateIn5x5; - for(size_t eleNr=0;eleNrsize();eleNr++){ auto elePtr = eleHandle->ptrAt(eleNr); eleTrkPtIso.push_back(calTrkIso(*elePtr,*eleHandle,candHandles)); @@ -230,6 +246,7 @@ void ElectronHEEPIDValueMapProducer::fillDescriptions(edm::ConfigurationDescript desc.add("eeRecHitsMiniAOD",edm::InputTag("reducedEcalRecHitsEE")); desc.add >("candsMiniAOD",{edm::InputTag("packedCandidates")}); desc.add("elesMiniAOD",edm::InputTag("gedGsfElectrons")); + desc.add("dataFormat",0); desc.add("trkIsoConfig",EleTkIsolFromCands::pSetDescript()); diff --git a/RecoEgamma/ElectronIdentification/python/heepIdVarValueMapProducer_cfi.py b/RecoEgamma/ElectronIdentification/python/heepIdVarValueMapProducer_cfi.py index 7d2928875a4ab..0e2e8a5c718d4 100644 --- a/RecoEgamma/ElectronIdentification/python/heepIdVarValueMapProducer_cfi.py +++ b/RecoEgamma/ElectronIdentification/python/heepIdVarValueMapProducer_cfi.py @@ -12,6 +12,7 @@ candsMiniAOD=cms.VInputTag("packedPFCandidates", "lostTracks"), elesMiniAOD=cms.InputTag("slimmedElectrons"), + dataFormat=cms.int32(0),#0 = auto detection, 1 = AOD, 2 = miniAOD trkIsoConfig= cms.PSet( barrelCuts=cms.PSet(