diff --git a/PhysicsTools/TagAndProbe/src/TagProbePairMaker.cc b/PhysicsTools/TagAndProbe/src/TagProbePairMaker.cc index c717a12d817c6..054dbb980fb16 100644 --- a/PhysicsTools/TagAndProbe/src/TagProbePairMaker.cc +++ b/PhysicsTools/TagAndProbe/src/TagProbePairMaker.cc @@ -120,6 +120,7 @@ tnp::TagProbePairMaker::arbitrate(TagProbePairs &pairs) const bool TTpair=false; for (TagProbePairs::iterator it2 = pairs.begin(); it2 != ed; ++it2) { // first check for Tag-Tag pairs + if (it2->tag.isNull()) continue; // skip already invalidated pairs if(it!=it2 && it->probe==it2->tag && it->tag==it2->probe){ //std::cout << "----------> probe is tag!!!!" << std::endl; TTpair=true; diff --git a/PhysicsTools/TagAndProbe/test/utilities/dumpPileup.py b/PhysicsTools/TagAndProbe/test/utilities/dumpPileup.py new file mode 100644 index 0000000000000..76b4771565c0a --- /dev/null +++ b/PhysicsTools/TagAndProbe/test/utilities/dumpPileup.py @@ -0,0 +1,7 @@ +#!/usr/bin/env python +import sys +import ROOT +fin=ROOT.TFile.Open(sys.argv[1]) +pu = fin.Get("pileup") +pileup = map( pu.GetBinContent, xrange(1,pu.GetNbinsX()+1) ) +print ",".join(map(lambda x: "%1.3g"%x, pileup)) diff --git a/RecoEgamma/ElectronIdentification/data/Fall17/effAreaElectrons_cone03_pfNeuHadronsAndPhotons_92X.txt b/RecoEgamma/ElectronIdentification/data/Fall17/effAreaElectrons_cone03_pfNeuHadronsAndPhotons_92X.txt new file mode 100644 index 0000000000000..709d8e2ca63ac --- /dev/null +++ b/RecoEgamma/ElectronIdentification/data/Fall17/effAreaElectrons_cone03_pfNeuHadronsAndPhotons_92X.txt @@ -0,0 +1,17 @@ +# This file contains Effective Area constants for +# computing pile-up corrections for the neutral hadron and photon +# isolation for an electron object. +# Documentation: +# +# https://indico.cern.ch/event/662749/contributions/2763091/attachments/1545124/2424854/talk_electron_ID_fall17.pdf +# +# The effective areas are based on 90% efficient contours +# +# |eta| min |eta| max effective area +0.0000 1.0000 0.1566 +1.0000 1.4790 0.1626 +1.4790 2.0000 0.1073 +2.0000 2.2000 0.0854 +2.2000 2.3000 0.1051 +2.3000 2.4000 0.1204 +2.4000 5.0000 0.1524 diff --git a/RecoEgamma/ElectronIdentification/plugins/cuts/GsfEleHadronicOverEMEnergyScaledCut.cc b/RecoEgamma/ElectronIdentification/plugins/cuts/GsfEleHadronicOverEMEnergyScaledCut.cc new file mode 100644 index 0000000000000..48df9132f8eee --- /dev/null +++ b/RecoEgamma/ElectronIdentification/plugins/cuts/GsfEleHadronicOverEMEnergyScaledCut.cc @@ -0,0 +1,64 @@ +#include "PhysicsTools/SelectorUtils/interface/CutApplicatorWithEventContentBase.h" +#include "DataFormats/EgammaCandidates/interface/GsfElectron.h" + +class GsfEleHadronicOverEMEnergyScaledCut : public CutApplicatorWithEventContentBase { +public: + GsfEleHadronicOverEMEnergyScaledCut(const edm::ParameterSet& c) : + CutApplicatorWithEventContentBase(c), + barrelC0_(c.getParameter("barrelC0")), + barrelCE_(c.getParameter("barrelCE")), + barrelCr_(c.getParameter("barrelCr")), + endcapC0_(c.getParameter("endcapC0")), + endcapCE_(c.getParameter("endcapCE")), + endcapCr_(c.getParameter("endcapCr")), + barrelCutOff_(c.getParameter("barrelCutOff")) + { + edm::InputTag rhoTag = c.getParameter("rho"); + contentTags_.emplace("rho",rhoTag); + + } + + result_type operator()(const reco::GsfElectronPtr&) const final; + + void setConsumes(edm::ConsumesCollector&) final; + void getEventContent(const edm::EventBase&) final; + + double value(const reco::CandidatePtr& cand) const final; + + CandidateType candidateType() const final { + return ELECTRON; + } + +private: + const float barrelC0_, barrelCE_, barrelCr_, endcapC0_, endcapCE_, endcapCr_, barrelCutOff_; + edm::Handle rhoHandle_; +}; + +DEFINE_EDM_PLUGIN(CutApplicatorFactory, + GsfEleHadronicOverEMEnergyScaledCut, + "GsfEleHadronicOverEMEnergyScaledCut"); + +void GsfEleHadronicOverEMEnergyScaledCut::setConsumes(edm::ConsumesCollector& cc) { + auto rho = cc.consumes(contentTags_["rho"]); + contentTokens_.emplace("rho",rho); +} + +void GsfEleHadronicOverEMEnergyScaledCut::getEventContent(const edm::EventBase& ev) { + ev.getByLabel(contentTags_["rho"],rhoHandle_); +} + + +CutApplicatorBase::result_type GsfEleHadronicOverEMEnergyScaledCut::operator()(const reco::GsfElectronPtr& cand) const { + + const double rho = rhoHandle_.isValid() ? (float)(*rhoHandle_) : 0; + const float energy = cand->superCluster()->energy(); + const float c0 = (std::abs(cand->superCluster()->position().eta()) < barrelCutOff_ ? barrelC0_ : endcapC0_); + const float cE = (std::abs(cand->superCluster()->position().eta()) < barrelCutOff_ ? barrelCE_ : endcapCE_); + const float cR = (std::abs(cand->superCluster()->position().eta()) < barrelCutOff_ ? barrelCr_ : endcapCr_); + return cand->hadronicOverEm() < c0 + cE/energy + cR*rho/energy; +} + +double GsfEleHadronicOverEMEnergyScaledCut::value(const reco::CandidatePtr& cand) const { + reco::GsfElectronPtr ele(cand); + return ele->hadronicOverEm(); +} diff --git a/RecoEgamma/ElectronIdentification/python/Identification/cutBasedElectronID_Fall17_94X_V1_cff.py b/RecoEgamma/ElectronIdentification/python/Identification/cutBasedElectronID_Fall17_94X_V1_cff.py new file mode 100644 index 0000000000000..a6bf3537616ce --- /dev/null +++ b/RecoEgamma/ElectronIdentification/python/Identification/cutBasedElectronID_Fall17_94X_V1_cff.py @@ -0,0 +1,190 @@ +from PhysicsTools.SelectorUtils.centralIDRegistry import central_id_registry + +import FWCore.ParameterSet.Config as cms + +# Common functions and classes for ID definition are imported here: +from RecoEgamma.ElectronIdentification.Identification.cutBasedElectronID_tools \ + import ( EleWorkingPoint_V4, + IsolationCutInputs_V2, + configureVIDCutBasedEleID_V4 ) + +# +# The ID cuts below are optimized IDs on Spring16 simulation with 80X-based production +# The cut values are taken from the twiki: +# https://twiki.cern.ch/twiki/bin/view/CMS/CutBasedElectronIdentificationRun2 +# (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): +# https://indico.cern.ch/event/662751/contributions/2778044/attachments/1562080/2459801/171121_egamma_workshop.pdf +# +# First, define cut values +# + +# Veto working point Barrel and Endcap +#V1 of IDs good for Moriond 18 +idName = "cutBasedElectronID-Fall17-94X-V1-veto" +WP_Veto_EB = EleWorkingPoint_V4( + idName = idName , # idName + full5x5_sigmaIEtaIEtaCut = 0.0128 , # full5x5_sigmaIEtaIEtaCut + dEtaInSeedCut = 0.00523 , # dEtaInSeedCut + dPhiInCut = 0.159 , # dPhiInCut + hOverECut_C0 = 0.05 , # hOverECut + hOverECut_CE = 1.12 , + hOverECut_Cr = 0.0368 , + relCombIsolationWithEALowPtCut = 0.168 , # relCombIsolationWithEALowPtCut + relCombIsolationWithEAHighPtCut= 0.168 , # relCombIsolationWithEAHighPtCut + absEInverseMinusPInverseCut = 0.193 , # absEInverseMinusPInverseCut + # conversion veto cut needs no parameters, so not mentioned + missingHitsCut = 2 # missingHitsCut + ) + +WP_Veto_EE = EleWorkingPoint_V4( + idName = idName , # idName + full5x5_sigmaIEtaIEtaCut = 0.0445 , # full5x5_sigmaIEtaIEtaCut + dEtaInSeedCut = 0.00984 , # dEtaInSeedCut + dPhiInCut = 0.157 , # dPhiInCut + hOverECut_C0 = 0.05 , # hOverECut + hOverECut_CE = 0.5 , + hOverECut_Cr = 0.201 , + relCombIsolationWithEALowPtCut = 0.185 , # relCombIsolationWithEALowPtCut + relCombIsolationWithEAHighPtCut= 0.185 , # relCombIsolationWithEAHighPtCut + absEInverseMinusPInverseCut = 0.0962 , # absEInverseMinusPInverseCut + # conversion veto cut needs no parameters, so not mentioned + missingHitsCut = 3 # missingHitsCut + ) + +# Loose working point Barrel and Endcap +idName = "cutBasedElectronID-Fall17-94X-V1-loose" +WP_Loose_EB = EleWorkingPoint_V4( + idName = idName , # idName + full5x5_sigmaIEtaIEtaCut = 0.0105 , # full5x5_sigmaIEtaIEtaCut + dEtaInSeedCut = 0.00387 , # dEtaInSeedCut + dPhiInCut = 0.0716 , # dPhiInCut + hOverECut_C0 = 0.05 , # hOverECut + hOverECut_CE = 1.12 , + hOverECut_Cr = 0.0368 , + relCombIsolationWithEALowPtCut = 0.133 , # relCombIsolationWithEALowPtCut + relCombIsolationWithEAHighPtCut= 0.133 , # relCombIsolationWithEAHighPtCut + absEInverseMinusPInverseCut = 0.129 , # absEInverseMinusPInverseCut + # conversion veto cut needs no parameters, so not mentioned + missingHitsCut = 1 # missingHitsCut + ) + +WP_Loose_EE = EleWorkingPoint_V4( + idName = idName , # idName + full5x5_sigmaIEtaIEtaCut = 0.0356 , # full5x5_sigmaIEtaIEtaCut + dEtaInSeedCut = 0.0072 , # dEtaInSeedCut + dPhiInCut = 0.147 , # dPhiInCut + hOverECut_C0 = 0.0414 , # hOverECut + hOverECut_CE = 0.5 , + hOverECut_Cr = 0.201 , + relCombIsolationWithEALowPtCut = 0.146 , # relCombIsolationWithEALowPtCut + relCombIsolationWithEAHighPtCut= 0.146 , # relCombIsolationWithEAHighPtCut + absEInverseMinusPInverseCut = 0.0875 , # absEInverseMinusPInverseCut + # conversion veto cut needs no parameters, so not mentioned + missingHitsCut = 1 # missingHitsCut + ) + +# Medium working point Barrel and Endcap +idName = "cutBasedElectronID-Fall17-94X-V1-medium" +WP_Medium_EB = EleWorkingPoint_V4( + idName = idName , # idName + full5x5_sigmaIEtaIEtaCut = 0.0105, # full5x5_sigmaIEtaIEtaCut + dEtaInSeedCut = 0.00365, # dEtaInSeedCut + dPhiInCut = 0.0588 , # dPhiInCut + hOverECut_C0 = 0.026 , # hOverECut + hOverECut_CE = 1.12 , + hOverECut_Cr = 0.0368 , + relCombIsolationWithEALowPtCut = 0.0718 , # relCombIsolationWithEALowPtCut + relCombIsolationWithEAHighPtCut= 0.0718 , # relCombIsolationWithEAHighPtCut + absEInverseMinusPInverseCut = 0.0327 , # absEInverseMinusPInverseCut + # conversion veto cut needs no parameters, so not mentioned + missingHitsCut = 1 # missingHitsCut + ) + +WP_Medium_EE = EleWorkingPoint_V4( + idName = idName , # idName + full5x5_sigmaIEtaIEtaCut = 0.0309 , # full5x5_sigmaIEtaIEtaCut + dEtaInSeedCut = 0.00625 , # dEtaInSeedCut + dPhiInCut = 0.0355 , # dPhiInCut + hOverECut_C0 = 0.026 , # hOverECut + hOverECut_CE = 0.5 , + hOverECut_Cr = 0.201 , + relCombIsolationWithEALowPtCut = 0.143 , # relCombIsolationWithEALowPtCut + relCombIsolationWithEAHighPtCut= 0.143 , # relCombIsolationWithEAHighPtCut + absEInverseMinusPInverseCut = 0.0335 , # absEInverseMinusPInverseCut + # conversion veto cut needs no parameters, so not mentioned + missingHitsCut = 1 # missingHitsCut + ) + +# Tight working point Barrel and Endcap +idName = "cutBasedElectronID-Fall17-94X-V1-tight" +WP_Tight_EB = EleWorkingPoint_V4( + idName = idName , # idName + full5x5_sigmaIEtaIEtaCut = 0.0104 , # full5x5_sigmaIEtaIEtaCut + dEtaInSeedCut = 0.00353 , # dEtaInSeedCut + dPhiInCut = 0.0499 , # dPhiInCut + hOverECut_C0 = 0.026 , # hOverECut + hOverECut_CE = 1.12 , + hOverECut_Cr = 0.0368 , + relCombIsolationWithEALowPtCut = 0.0361 , # relCombIsolationWithEALowPtCut + relCombIsolationWithEAHighPtCut= 0.0361 , # relCombIsolationWithEAHighPtCut + absEInverseMinusPInverseCut = 0.0278 , # absEInverseMinusPInverseCut + # conversion veto cut needs no parameters, so not mentioned + missingHitsCut = 1 # missingHitsCut + ) + +WP_Tight_EE = EleWorkingPoint_V4( + idName = idName , # idName + full5x5_sigmaIEtaIEtaCut = 0.0305 , # full5x5_sigmaIEtaIEtaCut + dEtaInSeedCut = 0.00567 , # dEtaInSeedCut + dPhiInCut = 0.0165 , # dPhiInCut + hOverECut_C0 = 0.026 , # hOverECut + hOverECut_CE = 0.5 , + hOverECut_Cr = 0.201 , + relCombIsolationWithEALowPtCut = 0.094 , # relCombIsolationWithEALowPtCut + relCombIsolationWithEAHighPtCut= 0.094 , # relCombIsolationWithEAHighPtCut + absEInverseMinusPInverseCut = 0.0158 , # absEInverseMinusPInverseCut + # conversion veto cut needs no parameters, so not mentioned + missingHitsCut = 1 # missingHitsCut + ) + +# Second, define what effective areas to use for pile-up correction +isoInputs = IsolationCutInputs_V2( + # phoIsolationEffAreas + "RecoEgamma/ElectronIdentification/data/Fall17/effAreaElectrons_cone03_pfNeuHadronsAndPhotons_92X.txt" +) + + +# +# Set up VID configuration for all cuts and working points +# + +cutBasedElectronID_Fall17_94X_V1_veto = configureVIDCutBasedEleID_V4(WP_Veto_EB, WP_Veto_EE, isoInputs) +cutBasedElectronID_Fall17_94X_V1_loose = configureVIDCutBasedEleID_V4(WP_Loose_EB, WP_Loose_EE, isoInputs) +cutBasedElectronID_Fall17_94X_V1_medium = configureVIDCutBasedEleID_V4(WP_Medium_EB, WP_Medium_EE, isoInputs) +cutBasedElectronID_Fall17_94X_V1_tight = configureVIDCutBasedEleID_V4(WP_Tight_EB, WP_Tight_EE, isoInputs) + + +# 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 "calculateIdMD5 +# 3) update the MD5 sum strings below and uncomment the lines again. +# + +central_id_registry.register(cutBasedElectronID_Fall17_94X_V1_veto.idName, + '43be9b381a8d9b0910b7f81a5ad8ff3a') +central_id_registry.register(cutBasedElectronID_Fall17_94X_V1_loose.idName, + '0b8456d622494441fe713a6858e0f7c1') +central_id_registry.register(cutBasedElectronID_Fall17_94X_V1_medium.idName, + 'a238ee70910de53d36866e89768500e9') +central_id_registry.register(cutBasedElectronID_Fall17_94X_V1_tight.idName, + '4acb2d2796efde7fba75380ce8823fc2') + + +### for now until we have a database... +cutBasedElectronID_Fall17_94X_V1_veto.isPOGApproved = cms.untracked.bool(True) +cutBasedElectronID_Fall17_94X_V1_loose.isPOGApproved = cms.untracked.bool(True) +cutBasedElectronID_Fall17_94X_V1_medium.isPOGApproved = cms.untracked.bool(True) +cutBasedElectronID_Fall17_94X_V1_tight.isPOGApproved = cms.untracked.bool(True) diff --git a/RecoEgamma/ElectronIdentification/python/Identification/cutBasedElectronID_tools.py b/RecoEgamma/ElectronIdentification/python/Identification/cutBasedElectronID_tools.py index f968e136a7371..70bae8fc1b676 100644 --- a/RecoEgamma/ElectronIdentification/python/Identification/cutBasedElectronID_tools.py +++ b/RecoEgamma/ElectronIdentification/python/Identification/cutBasedElectronID_tools.py @@ -112,6 +112,41 @@ def __init__(self, # conversion veto cut needs no parameters, so not mentioned self.missingHitsCut = missingHitsCut +class EleWorkingPoint_V4: + """ + This is a container class to hold numerical cut values for either + the barrel or endcap set of cuts for electron cut-based ID + With respect to V3, the hOverE cut is made energy and pileup dependend as presented in + https://indico.cern.ch/event/662749/contributions/2763092/attachments/1545209/2425054/talk_electron_ID_2017.pdf + """ + def __init__(self, + idName, + dEtaInSeedCut, + dPhiInCut, + full5x5_sigmaIEtaIEtaCut, + hOverECut_C0, + hOverECut_CE, + hOverECut_Cr, + absEInverseMinusPInverseCut, + relCombIsolationWithEALowPtCut, + relCombIsolationWithEAHighPtCut, + # conversion veto cut needs no parameters, so not mentioned + missingHitsCut + ): + self.idName = idName + self.dEtaInSeedCut = dEtaInSeedCut + self.dPhiInCut = dPhiInCut + self.full5x5_sigmaIEtaIEtaCut = full5x5_sigmaIEtaIEtaCut + self.hOverECut_C0 = hOverECut_C0 + self.hOverECut_CE = hOverECut_CE + self.hOverECut_Cr = hOverECut_Cr + self.absEInverseMinusPInverseCut = absEInverseMinusPInverseCut + self.relCombIsolationWithEALowPtCut = relCombIsolationWithEALowPtCut + self.relCombIsolationWithEAHighPtCut = relCombIsolationWithEAHighPtCut + # conversion veto cut needs no parameters, so not mentioned + self.missingHitsCut = missingHitsCut + + class EleHLTSelection_V1: """ This is a container class to hold numerical cut values for either @@ -225,6 +260,23 @@ def psetHadronicOverEMCut(wpEB, wpEE): isIgnored = cms.bool(False) ) +# Configure energy and pileup dependend H/E cut +def psetHadronicOverEMEnergyScaledCut(wpEB, wpEE): + return cms.PSet( + cutName = cms.string('GsfEleHadronicOverEMEnergyScaledCut'), + barrelC0 = cms.double( wpEB.hOverECut_C0 ), + barrelCE = cms.double( wpEB.hOverECut_CE ), + barrelCr = cms.double( wpEB.hOverECut_Cr ), + endcapC0 = cms.double( wpEE.hOverECut_C0 ), + endcapCE = cms.double( wpEE.hOverECut_CE ), + endcapCr = cms.double( wpEE.hOverECut_Cr ), + rho = cms.InputTag("fixedGridRhoFastjetAll"), + barrelCutOff = cms.double(ebCutOff), + needsAdditionalProducts = cms.bool(True), + isIgnored = cms.bool(False) + ) + + # Configure |1/E-1/p| cut def psetEInerseMinusPInverseCut(wpEB, wpEE): return cms.PSet( @@ -590,6 +642,35 @@ def configureVIDCutBasedEleID_V3( wpEB, wpEE, isoInputs ): # return parameterSet +def configureVIDCutBasedEleID_V4( wpEB, wpEE, isoInputs ): + """ + This function configures the full cms.PSet for a VID ID and returns it. + The inputs: two objects of the type WorkingPoint_V3, one + containing the cuts for the Barrel (EB) and the other one for the Endcap (EE). + The third argument is an object that contains information necessary + for isolation calculations. + In this version, the energy and pileup dependend hOverE is introduced + """ + # print "VID: Configuring cut set %s" % wpEB.idName + parameterSet = cms.PSet( + # + idName = cms.string( wpEB.idName ), # same name stored in the _EB and _EE objects + cutFlow = cms.VPSet( + psetMinPtCut(), + psetPhoSCEtaMultiRangeCut(), # eta cut + psetDEtaInSeedCut(wpEB, wpEE), # dEtaIn seed cut + psetDPhiInCut(wpEB, wpEE), # dPhiIn cut + psetPhoFull5x5SigmaIEtaIEtaCut(wpEB, wpEE), # full 5x5 sigmaIEtaIEta cut + psetHadronicOverEMEnergyScaledCut(wpEB, wpEE), # H/E cut + psetEInerseMinusPInverseCut(wpEB, wpEE), # |1/e-1/p| cut + psetEffAreaPFIsoCut(wpEB, wpEE, isoInputs), # rel. comb. PF isolation cut + psetConversionVetoCut(), + psetMissingHitsCut(wpEB, wpEE) + ) + ) + # + return parameterSet + # ----------------------------- # HLT-safe common definitions diff --git a/RecoEgamma/ElectronIdentification/test/runtests.sh b/RecoEgamma/ElectronIdentification/test/runtests.sh index 2648092eb8583..1850d07306351 100755 --- a/RecoEgamma/ElectronIdentification/test/runtests.sh +++ b/RecoEgamma/ElectronIdentification/test/runtests.sh @@ -1,14 +1,18 @@ function die { echo $1: status $2 ; exit $2; } ids_to_test=( -"RecoEgamma.ElectronIdentification.Identification.cutBasedElectronID_PHYS14_PU20bx25_V2_cff" -"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" + + "RecoEgamma.ElectronIdentification.Identification.cutBasedElectronID_PHYS14_PU20bx25_V2_cff" + "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" + 'RecoEgamma.ElectronIdentification.Identification.cutBasedElectronID_Summer16_80X_V1_cff' + 'RecoEgamma.ElectronIdentification.Identification.mvaElectronID_Spring16_GeneralPurpose_V1_cff' + 'RecoEgamma.ElectronIdentification.Identification.cutBasedElectronID_Fall17_94X_V1_cff' ) for id_set in "${ids_to_test[@]}"; do diff --git a/RecoEgamma/PhotonIdentification/data/Fall17/effAreaPhotons_cone03_pfChargedHadrons_90percentBased.txt b/RecoEgamma/PhotonIdentification/data/Fall17/effAreaPhotons_cone03_pfChargedHadrons_90percentBased.txt new file mode 100644 index 0000000000000..355d3fce259be --- /dev/null +++ b/RecoEgamma/PhotonIdentification/data/Fall17/effAreaPhotons_cone03_pfChargedHadrons_90percentBased.txt @@ -0,0 +1,16 @@ +# This file contains Effective Area constants for +# computing pile-up corrections for the charged hadron isolation +# for a photon object. +# +# The constants are based on 90% contours of isolation. +# These effective areas are also found on the following slides: +# https://indico.cern.ch/event/662751/contributions/2778043/attachments/1562017/2459674/EGamma_WorkShop_21.11.17_Debabrata.pdf +# +# |eta| min |eta| max effective area +0.0000 1.0000 0.057 +1.0000 1.4790 0.052 +1.4790 2.0000 0.045 +2.0000 2.2000 0.044 +2.2000 2.3000 0.043 +2.3000 2.4000 0.038 +2.4000 5.0000 0.034 diff --git a/RecoEgamma/PhotonIdentification/data/Fall17/effAreaPhotons_cone03_pfChargedHadrons_90percentBased_TrueVtx.txt b/RecoEgamma/PhotonIdentification/data/Fall17/effAreaPhotons_cone03_pfChargedHadrons_90percentBased_TrueVtx.txt new file mode 100644 index 0000000000000..2aad9998736f2 --- /dev/null +++ b/RecoEgamma/PhotonIdentification/data/Fall17/effAreaPhotons_cone03_pfChargedHadrons_90percentBased_TrueVtx.txt @@ -0,0 +1,16 @@ +# This file contains Effective Area constants for +# computing pile-up corrections for the charged hadron isolation +# for a photon object. +# +# The constants are based on 90% contours of isolation. +# These effective areas are also found on the following slides: +# https://indico.cern.ch/event/662751/contributions/2778043/attachments/1562017/2459674/EGamma_WorkShop_21.11.17_Debabrata.pdf +# +# |eta| min |eta| max effective area +0.0000 1.0000 0.038 +1.0000 1.4790 0.047 +1.4790 2.0000 0.043 +2.0000 2.2000 0.038 +2.2000 2.3000 0.033 +2.3000 2.4000 0.031 +2.4000 5.0000 0.027 diff --git a/RecoEgamma/PhotonIdentification/data/Fall17/effAreaPhotons_cone03_pfNeutralHadrons_90percentBased.txt b/RecoEgamma/PhotonIdentification/data/Fall17/effAreaPhotons_cone03_pfNeutralHadrons_90percentBased.txt new file mode 100644 index 0000000000000..aa37455b06aa3 --- /dev/null +++ b/RecoEgamma/PhotonIdentification/data/Fall17/effAreaPhotons_cone03_pfNeutralHadrons_90percentBased.txt @@ -0,0 +1,16 @@ +# This file contains Effective Area constants for +# computing pile-up corrections for the neutral hadron isolation +# for a photon object. +# +# The effective areas are derived based on the 90% contour method. +# These effective areas are also found on the following slides: +# https://indico.cern.ch/event/662751/contributions/2778043/attachments/1562017/2459674/EGamma_WorkShop_21.11.17_Debabrata.pdf +# +# |eta| min |eta| max effective area +0.0000 1.0000 0.062 +1.0000 1.4790 0.112 +1.4790 2.0000 0.077 +2.0000 2.2000 0.033 +2.2000 2.3000 0.007 +2.3000 2.4000 0.007 +2.4000 5.0000 0.012 diff --git a/RecoEgamma/PhotonIdentification/data/Fall17/effAreaPhotons_cone03_pfNeutralHadrons_90percentBased_TrueVtx.txt b/RecoEgamma/PhotonIdentification/data/Fall17/effAreaPhotons_cone03_pfNeutralHadrons_90percentBased_TrueVtx.txt new file mode 100644 index 0000000000000..117c377860841 --- /dev/null +++ b/RecoEgamma/PhotonIdentification/data/Fall17/effAreaPhotons_cone03_pfNeutralHadrons_90percentBased_TrueVtx.txt @@ -0,0 +1,16 @@ +# This file contains Effective Area constants for +# computing pile-up corrections for the neutral hadron isolation +# for a photon object. +# +# The effective areas are derived based on the 90% contour method. +# These effective areas are also found on the following slides: +# https://indico.cern.ch/event/662751/contributions/2778043/attachments/1562017/2459674/EGamma_WorkShop_21.11.17_Debabrata.pdf +# +# |eta| min |eta| max effective area +0.0000 1.0000 0.064 +1.0000 1.4790 0.110 +1.4790 2.0000 0.076 +2.0000 2.2000 0.023 +2.2000 2.3000 0.0151 +2.3000 2.4000 0.0001 +2.4000 5.0000 0.013 diff --git a/RecoEgamma/PhotonIdentification/data/Fall17/effAreaPhotons_cone03_pfPhotons_90percentBased.txt b/RecoEgamma/PhotonIdentification/data/Fall17/effAreaPhotons_cone03_pfPhotons_90percentBased.txt new file mode 100644 index 0000000000000..4068b0bcbe67a --- /dev/null +++ b/RecoEgamma/PhotonIdentification/data/Fall17/effAreaPhotons_cone03_pfPhotons_90percentBased.txt @@ -0,0 +1,16 @@ +# This file contains Effective Area constants for +# computing pile-up corrections for the photons isolation +# for a photon object. +# +# The constants are based on 90% contours of isolation. +# These effective areas are also found on the following slides: +# https://indico.cern.ch/event/662751/contributions/2778043/attachments/1562017/2459674/EGamma_WorkShop_21.11.17_Debabrata.pdf +# +# |eta| min |eta| max effective area +0.0000 1.0000 0.122 +1.0000 1.4790 0.110 +1.4790 2.0000 0.061 +2.0000 2.2000 0.077 +2.2000 2.3000 0.099 +2.3000 2.4000 0.109 +2.4000 5.0000 0.134 diff --git a/RecoEgamma/PhotonIdentification/data/Fall17/effAreaPhotons_cone03_pfPhotons_90percentBased_TrueVtx.txt b/RecoEgamma/PhotonIdentification/data/Fall17/effAreaPhotons_cone03_pfPhotons_90percentBased_TrueVtx.txt new file mode 100644 index 0000000000000..3ce289fd51c18 --- /dev/null +++ b/RecoEgamma/PhotonIdentification/data/Fall17/effAreaPhotons_cone03_pfPhotons_90percentBased_TrueVtx.txt @@ -0,0 +1,16 @@ +# This file contains Effective Area constants for +# computing pile-up corrections for the photons isolation +# for a photon object. +# +# The constants are based on 90% contours of isolation. +# These effective areas are also found on the following slides: +# https://indico.cern.ch/event/662751/contributions/2778043/attachments/1562017/2459674/EGamma_WorkShop_21.11.17_Debabrata.pdf +# +# |eta| min |eta| max effective area +0.0000 1.0000 0.124 +1.0000 1.4790 0.109 +1.4790 2.0000 0.063 +2.0000 2.2000 0.078 +2.2000 2.3000 0.100 +2.3000 2.4000 0.115 +2.4000 5.0000 0.137 diff --git a/RecoEgamma/PhotonIdentification/plugins/PhotonIDValueMapProducer.cc b/RecoEgamma/PhotonIdentification/plugins/PhotonIDValueMapProducer.cc index a17ad9cab6bd5..8aa8cdae8d2b1 100644 --- a/RecoEgamma/PhotonIdentification/plugins/PhotonIDValueMapProducer.cc +++ b/RecoEgamma/PhotonIdentification/plugins/PhotonIDValueMapProducer.cc @@ -67,7 +67,7 @@ class PhotonIDValueMapProducer : public edm::stream::EDProducer<> { float computeWorstPFChargedIsolation(const T& photon, const U& pfCandidates, const edm::Handle vertices, - bool isAOD, + bool isAOD, bool isPVConstraint,const reco::Vertex& pv, float dRmax, float dxyMax, float dzMax, float dRvetoBarrel, float dRvetoEndcap, float ptMin); @@ -76,7 +76,10 @@ class PhotonIDValueMapProducer : public edm::stream::EDProducer<> { reco::PFCandidate::ParticleType candidatePdgId(const edm::Ptr candidate, bool isAOD); - std::pair getTrackDxyDz(const edm::Ptr & candidate, const reco::Particle::Point & vtxpos, bool isAOD); + const reco::Track* getTrackPointer(const edm::Ptr candidate, bool isAOD); + void getImpactParameters(const edm::Ptr& candidate, + bool isAOD, const reco::Vertex& pv, float &dxy, float &dz); + // The object that will compute 5x5 quantities @@ -117,6 +120,13 @@ class PhotonIDValueMapProducer : public edm::stream::EDProducer<> { constexpr static char phoPhotonIsolation_[] = "phoPhotonIsolation"; constexpr static char phoWorstChargedIsolation_[] = "phoWorstChargedIsolation"; constexpr static char phoWorstChargedIsolationWithConeVeto_[] = "phoWorstChargedIsolationWithConeVeto"; + constexpr static char phoWorstChargedIsolationWithPVConstraint_[] = "phoWorstChargedIsolationWithPVConstraint"; + constexpr static char phoWorstChargedIsolationWithConeVetoWithPVConstraint_[] = "phoWorstChargedIsolationWithConeVetoWithPVConstraint"; + + //PFCluster Isolation + constexpr static char phoTrkIsolation_[] = "phoTrkIsolation"; + constexpr static char phoHcalPFClIsolation_[] = "phoHcalPFClIsolation"; + constexpr static char phoEcalPFClIsolation_[] = "phoEcalPFClIsolation"; }; @@ -134,6 +144,14 @@ constexpr char PhotonIDValueMapProducer::phoNeutralHadronIsolation_[]; constexpr char PhotonIDValueMapProducer::phoPhotonIsolation_[]; constexpr char PhotonIDValueMapProducer::phoWorstChargedIsolation_[]; constexpr char PhotonIDValueMapProducer::phoWorstChargedIsolationWithConeVeto_[]; +constexpr char PhotonIDValueMapProducer::phoWorstChargedIsolationWithPVConstraint_[]; +constexpr char PhotonIDValueMapProducer::phoWorstChargedIsolationWithConeVetoWithPVConstraint_[]; + +//PFCluster Isolation +constexpr char PhotonIDValueMapProducer::phoTrkIsolation_[]; +constexpr char PhotonIDValueMapProducer::phoHcalPFClIsolation_[]; +constexpr char PhotonIDValueMapProducer::phoEcalPFClIsolation_[]; + PhotonIDValueMapProducer::PhotonIDValueMapProducer(const edm::ParameterSet& iConfig) { @@ -197,6 +215,14 @@ PhotonIDValueMapProducer::PhotonIDValueMapProducer(const edm::ParameterSet& iCon produces >(phoPhotonIsolation_); produces >(phoWorstChargedIsolation_); produces >(phoWorstChargedIsolationWithConeVeto_); + produces >(phoWorstChargedIsolationWithPVConstraint_); + produces >(phoWorstChargedIsolationWithConeVetoWithPVConstraint_); + + //PFCluster Isolations + produces >(phoTrkIsolation_); + produces >(phoHcalPFClIsolation_); + produces >(phoEcalPFClIsolation_); + } @@ -292,6 +318,13 @@ void PhotonIDValueMapProducer::produce(edm::Event& iEvent, const edm::EventSetup std::vector phoPhotonIsolation; std::vector phoWorstChargedIsolation; std::vector phoWorstChargedIsolationWithConeVeto; + std::vector phoWorstChargedIsolationWithPVConstraint; + std::vector phoWorstChargedIsolationWithConeVetoWithPVConstraint; + + //PFCluster Isolations + std::vector phoTrkIsolation; + std::vector phoHcalPFClIsolation; + std::vector phoEcalPFClIsolation; // reco::Photon::superCluster() is virtual so we can exploit polymorphism for (unsigned idxpho = 0; idxpho < src->size(); ++idxpho) { @@ -329,6 +362,21 @@ void PhotonIDValueMapProducer::produce(edm::Event& iEvent, const edm::EventSetup iPho->superCluster()->y() - pv.y(), iPho->superCluster()->z() - pv.z()); + //PFCluster Isolations + phoTrkIsolation .push_back( iPho->trkSumPtSolidConeDR04()); + if (isAOD) + { + phoHcalPFClIsolation .push_back(0.f); + phoEcalPFClIsolation .push_back(0.f); + } + else + { + edm::Ptr patPhotonPtr(src->ptrAt(idxpho)); + phoHcalPFClIsolation .push_back(patPhotonPtr->hcalPFClusterIso()); + phoEcalPFClIsolation .push_back(patPhotonPtr->ecalPFClusterIso()); + } + + // Zero the isolation sums float chargedIsoSum = 0; float neutralHadronIsoSum = 0; @@ -378,10 +426,12 @@ void PhotonIDValueMapProducer::produce(edm::Event& iEvent, const edm::EventSetup if( thisCandidateType == reco::PFCandidate::h ){ // for charged hadrons, additionally check consistency // with the PV + float dxy = -999, dz=-999; + getImpactParameters(iCand, isAOD, pv, dxy, dz); + - auto dxydz = getTrackDxyDz(iCand, pv.position(), isAOD); - if ( fabs(dxydz.first) > dxyMax) continue; - if ( fabs(dxydz.second) > dzMax) continue; + if(fabs(dxy) > dxyMax) continue; + if (fabs(dz) > dzMax) continue; // The candidate is eligible, increment the isolaiton chargedIsoSum += iCand->pt(); @@ -403,9 +453,10 @@ void PhotonIDValueMapProducer::produce(edm::Event& iEvent, const edm::EventSetup float dRvetoBarrel = 0.0; float dRvetoEndcap = 0.0; float ptMin = 0.0; + bool isPVConstraint=false; float worstChargedIso = computeWorstPFChargedIsolation(iPho, pfCandidatesHandle, vertices, - isAOD, coneSizeDR, dxyMax, dzMax, + isAOD, isPVConstraint,pv,coneSizeDR, dxyMax, dzMax, dRvetoBarrel, dRvetoEndcap, ptMin); phoWorstChargedIsolation .push_back( worstChargedIso ); @@ -416,10 +467,28 @@ void PhotonIDValueMapProducer::produce(edm::Event& iEvent, const edm::EventSetup ptMin = 0.1; float worstChargedIsoWithConeVeto = computeWorstPFChargedIsolation(iPho, pfCandidatesHandle, vertices, - isAOD, coneSizeDR, dxyMax, dzMax, + isAOD, isPVConstraint,pv, coneSizeDR, dxyMax, dzMax, dRvetoBarrel, dRvetoEndcap, ptMin); phoWorstChargedIsolationWithConeVeto .push_back( worstChargedIsoWithConeVeto ); + isPVConstraint=true; + float worstChargedIsoWithPVConstraint = + computeWorstPFChargedIsolation(iPho, pfCandidatesHandle, vertices, + isAOD, isPVConstraint,pv,coneSizeDR, dxyMax, dzMax, + dRvetoBarrel, dRvetoEndcap, ptMin); + phoWorstChargedIsolationWithPVConstraint .push_back( worstChargedIsoWithPVConstraint ); + + // Worst isolation computed with cone vetos and a ptMin cut, as in + // Run 2 Hgg code. + dRvetoBarrel = 0.02; + dRvetoEndcap = 0.02; + ptMin = 0.1; + float worstChargedIsoWithConeVetoWithPVConstraint = + computeWorstPFChargedIsolation(iPho, pfCandidatesHandle, vertices, + isAOD,isPVConstraint, pv, coneSizeDR, dxyMax, dzMax, + dRvetoBarrel, dRvetoEndcap, ptMin); + phoWorstChargedIsolationWithConeVetoWithPVConstraint .push_back( worstChargedIsoWithConeVetoWithPVConstraint ); + } @@ -432,12 +501,19 @@ void PhotonIDValueMapProducer::produce(edm::Event& iEvent, const edm::EventSetup writeValueMap(iEvent, src, phoFull5x5E2x5Max, phoFull5x5E2x5Max_); writeValueMap(iEvent, src, phoFull5x5E5x5 , phoFull5x5E5x5_); writeValueMap(iEvent, src, phoESEffSigmaRR , phoESEffSigmaRR_); - // IsolationsOB + // Isolation writeValueMap(iEvent, src, phoChargedIsolation, phoChargedIsolation_); writeValueMap(iEvent, src, phoNeutralHadronIsolation, phoNeutralHadronIsolation_); writeValueMap(iEvent, src, phoPhotonIsolation, phoPhotonIsolation_); writeValueMap(iEvent, src, phoWorstChargedIsolation, phoWorstChargedIsolation_); writeValueMap(iEvent, src, phoWorstChargedIsolationWithConeVeto, phoWorstChargedIsolationWithConeVeto_); + writeValueMap(iEvent, src, phoWorstChargedIsolationWithPVConstraint, phoWorstChargedIsolationWithPVConstraint_); + writeValueMap(iEvent, src, phoWorstChargedIsolationWithConeVetoWithPVConstraint, phoWorstChargedIsolationWithConeVetoWithPVConstraint_); + //PFCluster Isolation + writeValueMap(iEvent, src, phoTrkIsolation, phoTrkIsolation_); + writeValueMap(iEvent, src, phoHcalPFClIsolation, phoHcalPFClIsolation_); + writeValueMap(iEvent, src, phoEcalPFClIsolation, phoEcalPFClIsolation_); + } void PhotonIDValueMapProducer::writeValueMap(edm::Event &iEvent, @@ -468,7 +544,7 @@ template float PhotonIDValueMapProducer ::computeWorstPFChargedIsolation(const T& photon, const U& pfCandidates, const edm::Handle vertices, - bool isAOD, + bool isAOD, bool isPVConstraint,const reco::Vertex& pv, float dRmax, float dxyMax, float dzMax, float dRvetoBarrel, float dRvetoEndcap, float ptMin){ @@ -504,10 +580,15 @@ ::computeWorstPFChargedIsolation(const T& photon, const U& pfCandidates, if (iCand->pt() < ptMin) continue; + + float dxy=-999, dz=-999; + if(isPVConstraint) getImpactParameters(iCand, isAOD, pv, dxy, dz); + else getImpactParameters(iCand, isAOD, *vtx, dxy, dz); + + - auto dxydz = getTrackDxyDz(iCand, vtx->position(), isAOD); - if ( fabs(dxydz.first) > dxyMax) continue; - if ( fabs(dxydz.second) > dzMax) continue; + if( fabs(dxy) > dxyMax) continue; + if ( fabs(dz) > dzMax) continue; float dR2 = deltaR2(photon_directionWrtVtx.Eta(), photon_directionWrtVtx.Phi(), iCand->eta(), iCand->phi()); @@ -525,6 +606,7 @@ ::computeWorstPFChargedIsolation(const T& photon, const U& pfCandidates, return worstIsolation; } + reco::PFCandidate::ParticleType PhotonIDValueMapProducer::candidatePdgId(const edm::Ptr candidate, bool isAOD){ @@ -546,15 +628,33 @@ PhotonIDValueMapProducer::candidatePdgId(const edm::Ptr candida return thisCandidateType; } -std::pair -PhotonIDValueMapProducer::getTrackDxyDz(const edm::Ptr & candidate, const reco::Particle::Point & vtxpos, bool isAOD) { +const reco::Track* +PhotonIDValueMapProducer::getTrackPointer(const edm::Ptr candidate, bool isAOD){ + + const reco::Track* theTrack = nullptr; + if( isAOD ) + theTrack = &*( ((const recoCandPtr) candidate)->trackRef()); + else + theTrack = &( ((const patCandPtr) candidate)->pseudoTrack()); + + return theTrack; +} + +void PhotonIDValueMapProducer::getImpactParameters(const edm::Ptr& candidate, + bool isAOD, const reco::Vertex& pv, + float &dxy, float &dz){ + dxy=-999; + dz=-999; if( isAOD ) { - const reco::Track & theTrack = *recoCandPtr(candidate)->trackRef(); - return std::make_pair(theTrack.dxy(vtxpos),theTrack.dz(vtxpos)); + const reco::Track *theTrack = &*( ((const recoCandPtr) candidate)->trackRef()); + dxy = theTrack->dxy(pv.position()); + dz = theTrack->dz(pv.position()); } else { - const pat::PackedCandidate & theCand = *(patCandPtr(candidate)); - return std::make_pair(theCand.dxy(vtxpos),theCand.dz(vtxpos)); + const pat::PackedCandidate & aCand = *(patCandPtr(candidate)); + dxy = aCand.dxy(pv.position()); + dz = aCand.dz(pv.position()); + } } diff --git a/RecoEgamma/PhotonIdentification/plugins/PhotonMVAEstimatorRunIIFall17.cc b/RecoEgamma/PhotonIdentification/plugins/PhotonMVAEstimatorRunIIFall17.cc new file mode 100644 index 0000000000000..33c5ab939f10d --- /dev/null +++ b/RecoEgamma/PhotonIdentification/plugins/PhotonMVAEstimatorRunIIFall17.cc @@ -0,0 +1,311 @@ +#include "RecoEgamma/PhotonIdentification/plugins/PhotonMVAEstimatorRunIIFall17.h" +#include "FWCore/ParameterSet/interface/FileInPath.h" +#include "TMVA/MethodBDT.h" +#include + +PhotonMVAEstimatorRunIIFall17::PhotonMVAEstimatorRunIIFall17(const edm::ParameterSet& conf): + AnyMVAEstimatorRun2Base(conf), + methodName_("BDTG method"), + phoChargedIsolationLabel_(conf.getParameter("phoChargedIsolation")), + phoPhotonIsolationLabel_(conf.getParameter("phoPhotonIsolation")), + phoWorstChargedIsolationLabel_(conf.getParameter("phoWorstChargedIsolation")), + rhoLabel_(conf.getParameter("rho")), + effectiveAreas_( (conf.getParameter("effAreasConfigFile")).fullPath()), + phoIsoPtScalingCoeff_(conf.getParameter>("phoIsoPtScalingCoeff")) + +{ + + // + // Construct the MVA estimators + // + tag_ = conf.getParameter("mvaTag"); + + const std::vector weightFileNames + = conf.getParameter >("weightFileNames"); + + if( weightFileNames.size() != nCategories ) + throw cms::Exception("MVA config failure: ") + << "wrong number of weightfiles" << std::endl; + + gbrForests_.clear(); + // The method name is just a key to retrieve this method later, it is not + // a control parameter for a reader (the full definition of the MVA type and + // everything else comes from the xml weight files). + + // Create a TMVA reader object for each category + for(uint i=0; i& particle, const edm::Event& iEvent) const { + + const int iCategory = findCategory( particle ); + + const std::vector vars = fillMVAVariables( particle, iEvent ) ; + const float result = gbrForests_.at(iCategory)->GetClassifier(vars.data()); + + // DEBUG + const bool debug = false; + // The list of variables here must match EXACTLY the list and order in the + // packMVAVariables() call for barrel and endcap in this file. + if( debug ){ + printf("Printout of the photon variable inputs for MVA:\n"); + printf(" varRawE_ %f\n", vars[0] ); + printf(" varR9_ %f\n", vars[1] ); + printf(" varSieie_ %f\n", vars[2] ); + printf(" varSCEtaWidth_ %f\n", vars[3] ); + printf(" varSCPhiWidth_ %f\n", vars[4] ); + printf(" varSieip_ %f\n", vars[5] ); + printf(" varE2x2overE5x5_ %f\n", vars[6] ); + printf(" varPhoIsoRaw_ %f\n", vars[7] ); + printf(" varChIsoRaw_ %f\n", vars[8] ); + printf(" varWorstChRaw_ %f\n", vars[9] ); + printf(" varSCEta_ %f\n", vars[10] ); + printf(" varRho_ %f\n", vars[11] ); + if( isEndcapCategory( iCategory ) ) { + printf(" varESEffSigmaRR_ %f\n", vars[12] ); // for endcap MVA only + printf(" varESEnOverRawE_ %f\n", vars[13] ); // for endcap MVA only + } + } + + return result; +} + +int PhotonMVAEstimatorRunIIFall17::findCategory( const edm::Ptr& particle) const { + + // Try to cast the particle into a reco particle. + // This should work for both reco and pat. + const edm::Ptr phoRecoPtr = ( edm::Ptr )particle; + if( phoRecoPtr.isNull() ) + throw cms::Exception("MVA failure: ") + << " given particle is expected to be reco::Photon or pat::Photon," << std::endl + << " but appears to be neither" << std::endl; + + float eta = phoRecoPtr->superCluster()->eta(); + + // + // Determine the category + // + int iCategory = UNDEFINED; + const float ebeeSplit = 1.479; // division between barrel and endcap + + if ( std::abs(eta) < ebeeSplit) + iCategory = CAT_EB; + else + iCategory = CAT_EE; + + return iCategory; +} + +bool PhotonMVAEstimatorRunIIFall17:: +isEndcapCategory(int category) const { + + // For this specific MVA the function is trivial, but kept for possible + // future evolution to an MVA with more categories in eta + bool isEndcap = false; + if( category == CAT_EE ) + isEndcap = true; + + return isEndcap; +} + + +std::unique_ptr PhotonMVAEstimatorRunIIFall17:: +createSingleReader(const int iCategory, const edm::FileInPath &weightFile){ + + // + // Create the reader + // + TMVA::Reader tmpTMVAReader( "!Color:Silent:Error" ); + + // + // Configure all variables and spectators. Note: the order and names + // must match what is found in the xml weights file! + // + + tmpTMVAReader.AddVariable("SCRawE", &allMVAVars_.varRawE); + tmpTMVAReader.AddVariable("r9", &allMVAVars_.varR9); + tmpTMVAReader.AddVariable("sigmaIetaIeta", &allMVAVars_.varSieie); + tmpTMVAReader.AddVariable("etaWidth", &allMVAVars_.varSCEtaWidth); + tmpTMVAReader.AddVariable("phiWidth", &allMVAVars_.varSCPhiWidth); + tmpTMVAReader.AddVariable("covIEtaIPhi", &allMVAVars_.varSieip); + tmpTMVAReader.AddVariable("s4", &allMVAVars_.varE2x2overE5x5); + tmpTMVAReader.AddVariable("phoIso03", &allMVAVars_.varPhoIsoRaw); + tmpTMVAReader.AddVariable("chgIsoWrtChosenVtx", &allMVAVars_.varChIsoRaw); + tmpTMVAReader.AddVariable("chgIsoWrtWorstVtx", &allMVAVars_.varWorstChRaw); + tmpTMVAReader.AddVariable("scEta", &allMVAVars_.varSCEta); + tmpTMVAReader.AddVariable("rho", &allMVAVars_.varRho); + + // Endcap only variables + if( isEndcapCategory(iCategory) ){ + tmpTMVAReader.AddVariable("esEffSigmaRR", &allMVAVars_.varESEffSigmaRR); + tmpTMVAReader.AddVariable("esEnergy/SCRawE", &allMVAVars_.varESEnOverRawE); + } + + // + // Book the method and set up the weights file + // + tmpTMVAReader.BookMVA(methodName_ , weightFile.fullPath() ); + + return std::unique_ptr( new GBRForest( dynamic_cast( tmpTMVAReader.FindMVA(methodName_) ) ) ); + +} + +// A function that should work on both pat and reco objects +std::vector PhotonMVAEstimatorRunIIFall17::fillMVAVariables(const edm::Ptr& particle, const edm::Event& iEvent) const { + + // + // Declare all value maps corresponding to the above tokens + // + edm::Handle > full5x5SigmaIEtaIEtaMap; + edm::Handle > full5x5SigmaIEtaIPhiMap; + edm::Handle > full5x5E2x2Map; + edm::Handle > full5x5E5x5Map; + edm::Handle > esEffSigmaRRMap; + // + edm::Handle > phoChargedIsolationMap; + edm::Handle > phoPhotonIsolationMap; + edm::Handle > phoWorstChargedIsolationMap; + + // Rho will be pulled from the event content + edm::Handle rho; + + // Get the isolation maps + iEvent.getByLabel(phoChargedIsolationLabel_, phoChargedIsolationMap); + iEvent.getByLabel(phoPhotonIsolationLabel_, phoPhotonIsolationMap); + iEvent.getByLabel(phoWorstChargedIsolationLabel_, phoWorstChargedIsolationMap); + + // Get rho + iEvent.getByLabel(rhoLabel_,rho); + + // Make sure everything is retrieved successfully + if(! ( phoChargedIsolationMap.isValid() + && phoPhotonIsolationMap.isValid() + && phoWorstChargedIsolationMap.isValid() + && rho.isValid() ) ) + throw cms::Exception("MVA failure: ") + << "Failed to retrieve event content needed for this MVA" + << std::endl + << "Check python MVA configuration file and make sure all needed" + << std::endl + << "producers are running upstream" << std::endl; + + // Try to cast the particle into a reco particle. + // This should work for both reco and pat. + const edm::Ptr phoRecoPtr(particle); + if( phoRecoPtr.isNull() ) + throw cms::Exception("MVA failure: ") + << " given particle is expected to be reco::Photon or pat::Photon," << std::endl + << " but appears to be neither" << std::endl; + + // Both pat and reco particles have exactly the same accessors. + auto superCluster = phoRecoPtr->superCluster(); + // Full 5x5 cluster shapes. We could take some of this directly from + // the photon object, but some of these are not available. + float e2x2 = std::numeric_limits::max(); + float e5x5 = std::numeric_limits::max(); + float full5x5_sigmaIetaIeta = std::numeric_limits::max(); + float full5x5_sigmaIetaIphi = std::numeric_limits::max(); + float effSigmaRR = std::numeric_limits::max(); + + AllVariables allMVAVars; + + const auto& full5x5_pss = phoRecoPtr->full5x5_showerShapeVariables(); + e2x2 = full5x5_pss.e2x2; + e5x5 = full5x5_pss.e5x5; + full5x5_sigmaIetaIeta = full5x5_pss.sigmaIetaIeta; + full5x5_sigmaIetaIphi = full5x5_pss.sigmaIetaIphi; + effSigmaRR = full5x5_pss.effSigmaRR; + + allMVAVars.varRawE = superCluster->rawEnergy(); + allMVAVars.varR9 = phoRecoPtr->r9() ; + allMVAVars.varSieie = full5x5_sigmaIetaIeta; + allMVAVars.varSCEtaWidth = superCluster->etaWidth(); + allMVAVars.varSCPhiWidth = superCluster->phiWidth(); + allMVAVars.varSieip = full5x5_sigmaIetaIphi; + allMVAVars.varE2x2overE5x5 = e2x2/e5x5; + allMVAVars.varPhoIsoRaw = (*phoPhotonIsolationMap)[phoRecoPtr]; + allMVAVars.varChIsoRaw = (*phoChargedIsolationMap)[phoRecoPtr]; + allMVAVars.varWorstChRaw = (*phoWorstChargedIsolationMap)[phoRecoPtr]; + allMVAVars.varSCEta = superCluster->eta(); + allMVAVars.varRho = *rho; + allMVAVars.varESEffSigmaRR = effSigmaRR; + allMVAVars.varESEnOverRawE = superCluster->preshowerEnergy() / superCluster->rawEnergy(); + + constrainMVAVariables(allMVAVars); + // + // Important: the order of variables in the "vars" vector has to be EXACTLY + // the same as in the .xml file defining the MVA. + // + std::vector vars; + if( isEndcapCategory( findCategory( particle ) ) ) { + vars = packMVAVariables( + allMVAVars.varRawE, + allMVAVars.varR9, + allMVAVars.varSieie, + allMVAVars.varSCEtaWidth, + allMVAVars.varSCPhiWidth, + allMVAVars.varSieip, + allMVAVars.varE2x2overE5x5, + allMVAVars.varPhoIsoRaw, + allMVAVars.varChIsoRaw, + allMVAVars.varWorstChRaw, + allMVAVars.varSCEta, + allMVAVars.varRho, + allMVAVars.varESEffSigmaRR, + allMVAVars.varESEnOverRawE + ) + ; + } else { + vars = packMVAVariables( + allMVAVars.varRawE, + allMVAVars.varR9, + allMVAVars.varSieie, + allMVAVars.varSCEtaWidth, + allMVAVars.varSCPhiWidth, + allMVAVars.varSieip, + allMVAVars.varE2x2overE5x5, + allMVAVars.varPhoIsoRaw, + allMVAVars.varChIsoRaw, + allMVAVars.varWorstChRaw, + allMVAVars.varSCEta, + allMVAVars.varRho + ) + ; + } + return vars; +} + +void PhotonMVAEstimatorRunIIFall17::constrainMVAVariables(AllVariables&)const { + + // Check that variables do not have crazy values + + // This function is currently empty as this specific MVA was not + // developed with restricting variables to specific physical ranges. + +} + +void PhotonMVAEstimatorRunIIFall17::setConsumes(edm::ConsumesCollector&& cc) const { + cc.consumes >(phoChargedIsolationLabel_); + cc.consumes >(phoPhotonIsolationLabel_); + cc.consumes >( phoWorstChargedIsolationLabel_); + cc.consumes(rhoLabel_); +} + +DEFINE_EDM_PLUGIN(AnyMVAEstimatorRun2Factory, + PhotonMVAEstimatorRunIIFall17, + "PhotonMVAEstimatorRunIIFall17"); + diff --git a/RecoEgamma/PhotonIdentification/plugins/PhotonMVAEstimatorRunIIFall17.h b/RecoEgamma/PhotonIdentification/plugins/PhotonMVAEstimatorRunIIFall17.h new file mode 100644 index 0000000000000..b4d34420bc793 --- /dev/null +++ b/RecoEgamma/PhotonIdentification/plugins/PhotonMVAEstimatorRunIIFall17.h @@ -0,0 +1,127 @@ +#ifndef RecoEgamma_PhotonIdentification_PhotonMVAEstimatorRunIIFall17_H +#define RecoEgamma_PhotonIdentification_PhotonMVAEstimatorRunIIFall17_H + +#include "FWCore/Framework/interface/ConsumesCollector.h" + +#include "DataFormats/Common/interface/ValueMap.h" + +#include "RecoEgamma/EgammaTools/interface/AnyMVAEstimatorRun2Base.h" + +#include "DataFormats/EgammaCandidates/interface/Photon.h" + +#include "CondFormats/EgammaObjects/interface/GBRForest.h" + +#include "RecoEgamma/EgammaTools/interface/EffectiveAreas.h" + +#include +#include + +#include "TMVA/Factory.h" +#include "TMVA/Tools.h" +#include "TMVA/Reader.h" +#include "TMVA/DataLoader.h" + +class PhotonMVAEstimatorRunIIFall17 : public AnyMVAEstimatorRun2Base{ + + public: + + // Define here the number and the meaning of the categories + // for this specific MVA + const uint nCategories = 2; + enum mvaCategories { + UNDEFINED = -1, + CAT_EB = 0, + CAT_EE = 1 + }; + + // Define the struct that contains all necessary for MVA variables + struct AllVariables { + + float varRawE; + float varR9; + float varSieie; + float varSCEtaWidth; + float varSCPhiWidth; + float varSieip; + float varE2x2overE5x5; + float varSCEta; + float varESEnOverRawE; // for endcap MVA only + float varESEffSigmaRR; // for endcap MVA only + // Pile-up + float varRho; + // Isolations + float varPhoIsoRaw; + float varChIsoRaw; + float varWorstChRaw; + + }; + + // Constructor and destructor + PhotonMVAEstimatorRunIIFall17(const edm::ParameterSet& conf); + ~PhotonMVAEstimatorRunIIFall17(); + + // Calculation of the MVA value + float mvaValue( const edm::Ptr& particle, const edm::Event&) const; + + // Utility functions + std::unique_ptr createSingleReader(const int iCategory, const edm::FileInPath &weightFile); + + virtual int getNCategories() const { return nCategories; } + bool isEndcapCategory( int category ) const; + virtual const std::string& getName() const override final { return name_; } + virtual const std::string& getTag() const override final { return tag_; } + + // Functions that should work on both pat and reco electrons + // (use the fact that pat::Electron inherits from reco::GsfElectron) + std::vector fillMVAVariables(const edm::Ptr& particle, const edm::Event& iEvent) const override; + int findCategory( const edm::Ptr& particle ) const override; + // The function below ensures that the variables passed to MVA are + // within reasonable bounds + void constrainMVAVariables(AllVariables&) const; + + // Call this function once after the constructor to declare + // the needed event content pieces to the framework + void setConsumes(edm::ConsumesCollector&&) const override; + // Call this function once per event to retrieve all needed + // event content pices + + + + private: + + // MVA name. This is a unique name for this MVA implementation. + // It will be used as part of ValueMap names. + // For simplicity, keep it set to the class name. + const std::string name_ = "PhotonMVAEstimatorRunIIFall17"; + + // MVA tag. This is an additional string variable to distinguish + // instances of the estimator of this class configured with different + // weight files. + std::string tag_; + + // Data members + std::vector< std::unique_ptr > gbrForests_; + + // All variables needed by this MVA + const std::string methodName_; + AllVariables allMVAVars_; + + // This MVA implementation relies on several ValueMap objects + // produced upstream. + + // + // Declare all tokens that will be needed to retrieve misc + // data from the event content required by this MVA + // + const edm::InputTag phoChargedIsolationLabel_; + const edm::InputTag phoPhotonIsolationLabel_; + const edm::InputTag phoWorstChargedIsolationLabel_; + const edm::InputTag rhoLabel_; + + // Other objects needed by the MVA + EffectiveAreas effectiveAreas_; + std::vector phoIsoPtScalingCoeff_; + +}; + +#endif diff --git a/RecoEgamma/PhotonIdentification/python/Identification/cutBasedPhotonID_Fall17_94X_V1_TrueVtx_cff.py b/RecoEgamma/PhotonIdentification/python/Identification/cutBasedPhotonID_Fall17_94X_V1_TrueVtx_cff.py new file mode 100644 index 0000000000000..04020906e273c --- /dev/null +++ b/RecoEgamma/PhotonIdentification/python/Identification/cutBasedPhotonID_Fall17_94X_V1_TrueVtx_cff.py @@ -0,0 +1,156 @@ + +from PhysicsTools.SelectorUtils.centralIDRegistry import central_id_registry + +import FWCore.ParameterSet.Config as cms + +# Common functions and classes for ID definition are imported here: +from RecoEgamma.PhotonIdentification.Identification.cutBasedPhotonID_tools \ + import ( WorkingPoint_V2, + IsolationCutInputs, + configureVIDCutBasedPhoID_V5 ) + +# +# This is the first version of Spring16 cuts for 80X samples +# +# The cut values are taken from the twiki: +# https://twiki.cern.ch/twiki/bin/viewauth/CMS/CutBasedPhotonIdentificationRun2 +# (where they may not stay, if a newer version of cuts becomes available for these +# conditions) +# See also the presentation explaining these working points : +# https://indico.cern.ch/event/662751/contributions/2778043/attachments/1562017/2459674/EGamma_WorkShop_21.11.17_Debabrata.pdf + +# +# First, define cut values +# + +# Loose working point Barrel and Endcap +idName = "cutBasedPhotonID-Fall17-94X-V1-loose" +WP_Loose_EB = WorkingPoint_V2( + idName , # idName + 0.105 , # hOverECut + 0.0103 , # full5x5_SigmaIEtaIEtaCut +# Isolation cuts are generally absIso < C1 + pt*C2, except for NeuHad is < C1 + pt*C2 + pt*pt*C3 + 2.839 , # absPFChaHadIsoWithEACut_C1 + 0 , # absPFChaHadIsoWithEACut_C2 + 9.188 , # absPFNeuHadIsoWithEACut_C1 + 0.0126 , # absPFNeuHadIsoWithEACut_C2 #################check for AllPV + 0.000026 , # absPFNeuHadIsoWithEACut_C3 + 2.956 , # absPFPhoIsoWithEACut_C1 + 0.0035 # absPFPhoIsoWithEACut_C2 + ) +WP_Loose_EE = WorkingPoint_V2( + idName , #idName + 0.029 , # hOverECut + 0.0276 , # full5x5_SigmaIEtaIEtaCut +# Isolation cuts are generally absIso < C1 + pt*C2, except for NeuHad is < C1 + pt*C2 + pt*pt*C3 + 2.150 , # absPFChaHadIsoWithEACut_C1 + 0.00 , # absPFChaHadIsoWithEACut_C2 + 10.471 , # absPFNeuHadIsoWithEACut_C1 + 0.0119 , # absPFNeuHadIsoWithEACut_C2 + 0.000025 , # absPFNeuHadIsoWithEACut_C3 + 4.895 , # absPFPhoIsoWithEACut_C1 + 0.0040 # absPFPhoIsoWithEACut_C2 + ) + +# Medium working point Barrel and Endcap +idName = "cutBasedPhotonID-Fall17-94X-V1-medium" +WP_Medium_EB = WorkingPoint_V2( + idName , # idName + 0.035 , # hOverECut + 0.0103 , # full5x5_SigmaIEtaIEtaCut +# Isolation cuts are generally absIso < C1 + pt*C2, except for NeuHad is < C1 + pt*C2 + pt*pt*C3 + 1.416 , # absPFChaHadIsoWithEACut_C1 + 0.00 , # absPFChaHadIsoWithEACut_C2 + 2.491 , # absPFNeuHadIsoWithEACut_C1 + 0.0126 , # absPFNeuHadIsoWithEACut_C2 + 0.000026 , # absPFNeuHadIsowithEACut_C3 + 2.952 , # absPFPhoIsoWithEACut_C1 + 0.0040 # absPFPhoIsoWithEACut_C2 + ) + +WP_Medium_EE = WorkingPoint_V2( #################check for AllPV + idName , #idName + 0.027 , # hOverECut + 0.0271 , # full5x5_SigmaIEtaIEtaCut +# Isolation cuts are generally absIso < C1 + pt*C2, except for NeuHad is < C1 + pt*C2 + pt*pt*C3 + 1.012 , # absPFChaHadIsoWithEACut_C1 + 0.00 , # absPFChaHadIsoWithEACut_C2 + 9.131 , # absPFNeuHadIsoWithEACut_C1 + 0.0119 , # absPFNeuHadIsoWithEACut_C2 + 0.000025 , # absPFNeuHadIsowithEACut_C3 + 4.095 , # absPFPhoIsoWithEACut_C1 + 0.0040 # absPFPhoIsoWithEACut_C2 + ) + +# Tight working point Barrel and Endcap +idName = "cutBasedPhotonID-Fall17-94X-V1-tight" +WP_Tight_EB = WorkingPoint_V2( + idName , # idName + 0.020 , # hOverECut + 0.0103 , # full5x5_SigmaIEtaIEtaCut +# Isolation cuts are generally absIso < C1 + pt*C2, except for NeuHad is < C1 + pt*C2 + pt*pt*C3 + 1.158 , # absPFChaHadIsoWithEACut_C1 + 0.00 , # absPFChaHadIsoWithEACut_C2 + 1.267 , # absPFNeuHadIsoWithEACut_C1 + 0.0126 , # absPFNeuHadIsoWithEACut_C2 + 0.000026 , # absPFNeuHadIsoWithEACut_C3 + 2.065 , # absPFPhoIsoWithEACut_C1 + 0.0035 # absPFPhoIsoWithEACut_C2 + ) + +WP_Tight_EE = WorkingPoint_V2( + idName , #idName + 0.025 , # hOverECut + 0.0271 , # full5x5_SigmaIEtaIEtaCut +# Isolation cuts are generally absIso < C1 + pt*C2, except for NeuHad is < C1 + pt*C2 + pt*pt*C3 + 0.575 , # absPFChaHadIsoWithEACut_C1 + 0.00 , # absPFChaHadIsoWithEACut_C2 + 8.916 , # absPFNeuHadIsoWithEACut_C1 + 0.0119 , # absPFNeuHadIsoWithEACut_C2 + 0.000025 , # absPFNeuHadIsowithEACut_C3 + 3.272 , # absPFPhoIsoWithEACut_C1 + 0.0040 # absPFPhoIsoWithEACut_C2 + ) + + +# Second, define where to find the precomputed isolations and what effective +# areas to use for pile-up correction +isoInputs = IsolationCutInputs( + # chHadIsolationMapName + 'photonIDValueMapProducer:phoChargedIsolation' , + # chHadIsolationEffAreas + "RecoEgamma/PhotonIdentification/data/Fall17/effAreaPhotons_cone03_pfChargedHadrons_90percentBased_TrueVtx.txt", + # neuHadIsolationMapName + 'photonIDValueMapProducer:phoNeutralHadronIsolation' , + # neuHadIsolationEffAreas + "RecoEgamma/PhotonIdentification/data/Fall17/effAreaPhotons_cone03_pfNeutralHadrons_90percentBased_TrueVtx.txt" , + # phoIsolationMapName + "photonIDValueMapProducer:phoPhotonIsolation" , + # phoIsolationEffAreas + "RecoEgamma/PhotonIdentification/data/Fall17/effAreaPhotons_cone03_pfPhotons_90percentBased_TrueVtx.txt" +) + +# +# Finally, set up VID configuration for all cuts +# +cutBasedPhotonID_Fall17_94X_V1_loose = configureVIDCutBasedPhoID_V5 ( WP_Loose_EB, WP_Loose_EE, isoInputs) +cutBasedPhotonID_Fall17_94X_V1_medium = configureVIDCutBasedPhoID_V5 ( WP_Medium_EB, WP_Medium_EE, isoInputs) +cutBasedPhotonID_Fall17_94X_V1_tight = configureVIDCutBasedPhoID_V5 ( WP_Tight_EB, WP_Tight_EE, isoInputs) + +## 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(cutBasedPhotonID_Fall17_94X_V1_loose.idName, + '45515ee95e01fa36972ff7ba69186c97') +central_id_registry.register(cutBasedPhotonID_Fall17_94X_V1_medium.idName, + '772f7921fa146b630e4dbe79e475a421') +central_id_registry.register(cutBasedPhotonID_Fall17_94X_V1_tight.idName, + 'e260fee6f9011fb13ff56d45cccd21c5') + +cutBasedPhotonID_Fall17_94X_V1_loose.isPOGApproved = cms.untracked.bool(True) +cutBasedPhotonID_Fall17_94X_V1_medium.isPOGApproved = cms.untracked.bool(True) +cutBasedPhotonID_Fall17_94X_V1_tight.isPOGApproved = cms.untracked.bool(True) diff --git a/RecoEgamma/PhotonIdentification/python/Identification/cutBasedPhotonID_Fall17_94X_V1_cff.py b/RecoEgamma/PhotonIdentification/python/Identification/cutBasedPhotonID_Fall17_94X_V1_cff.py new file mode 100644 index 0000000000000..d18dff7ffdd79 --- /dev/null +++ b/RecoEgamma/PhotonIdentification/python/Identification/cutBasedPhotonID_Fall17_94X_V1_cff.py @@ -0,0 +1,155 @@ + +from PhysicsTools.SelectorUtils.centralIDRegistry import central_id_registry + +import FWCore.ParameterSet.Config as cms + +# Common functions and classes for ID definition are imported here: +from RecoEgamma.PhotonIdentification.Identification.cutBasedPhotonID_tools \ + import ( WorkingPoint_V2, + IsolationCutInputs, + configureVIDCutBasedPhoID_V5 ) + +# +# This is the first version of Spring16 cuts for 80X samples +# +# The cut values are taken from the twiki: +# https://twiki.cern.ch/twiki/bin/viewauth/CMS/CutBasedPhotonIdentificationRun2 +# (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): +# https://indico.cern.ch/event/662751/contributions/2778043/attachments/1562017/2459674/EGamma_WorkShop_21.11.17_Debabrata.pdf + +# +# First, define cut values +# + +# Loose working point Barrel and Endcap +idName = "cutBasedPhotonID-Fall17-94X-V1-loose" +WP_Loose_EB = WorkingPoint_V2( + idName , # idName + 0.043 , # hOverECut + 0.0101 , # full5x5_SigmaIEtaIEtaCut +# Isolation cuts are generally absIso < C1 + pt*C2, except for NeuHad is < C1 + pt*C2 + pt*pt*C3 + 1.403 , # absPFChaHadIsoWithEACut_C1 + 0 , # absPFChaHadIsoWithEACut_C2 + 15.959 , # absPFNeuHadIsoWithEACut_C1 + 0.0127 , # absPFNeuHadIsoWithEACut_C2 + 0.000026 , # absPFNeuHadIsoWithEACut_C3 + 3.06 , # absPFPhoIsoWithEACut_C1 + 0.0038 # absPFPhoIsoWithEACut_C2 + ) +WP_Loose_EE = WorkingPoint_V2( + idName , #idName + 0.026 , # hOverECut + 0.0267 , # full5x5_SigmaIEtaIEtaCut +# Isolation cuts are generally absIso < C1 + pt*C2, except for NeuHad is < C1 + pt*C2 + pt*pt*C3 + 2.809 , # absPFChaHadIsoWithEACut_C1 + 0.00 , # absPFChaHadIsoWithEACut_C2 + 7.056 , # absPFNeuHadIsoWithEACut_C1 + 0.0117 , # absPFNeuHadIsoWithEACut_C2 + 0.000025 , # absPFNeuHadIsoWithEACut_C3 + 4.766 , # absPFPhoIsoWithEACut_C1 + 0.0038 # absPFPhoIsoWithEACut_C2 + ) + +# Medium working point Barrel and Endcap +idName = "cutBasedPhotonID-Fall17-94X-V1-medium" +WP_Medium_EB = WorkingPoint_V2( + idName , # idName + 0.032 , # hOverECut + 0.0101 , # full5x5_SigmaIEtaIEtaCut +# Isolation cuts are generally absIso < C1 + pt*C2, except for NeuHad is < C1 + pt*C2 + pt*pt*C3 + 0.43 , # absPFChaHadIsoWithEACut_C1 + 0.00 , # absPFChaHadIsoWithEACut_C2 + 2.133 , # absPFNeuHadIsoWithEACut_C1 + 0.0127 , # absPFNeuHadIsoWithEACut_C2 + 0.000026 , # absPFNeuHadIsowithEACut_C3 + 2.344 , # absPFPhoIsoWithEACut_C1 + 0.0038 # absPFPhoIsoWithEACut_C2 + ) + +WP_Medium_EE = WorkingPoint_V2( + idName , #idName + 0.0219 , # hOverECut + 0.03001 , # full5x5_SigmaIEtaIEtaCut +# Isolation cuts are generally absIso < C1 + pt*C2, except for NeuHad is < C1 + pt*C2 + pt*pt*C3 + 0.442 , # absPFChaHadIsoWithEACut_C1 + 0.00 , # absPFChaHadIsoWithEACut_C2 + 1.715 , # absPFNeuHadIsoWithEACut_C1 + 0.0117 , # absPFNeuHadIsoWithEACut_C2 + 0.000025 , # absPFNeuHadIsowithEACut_C3 + 3.863 , # absPFPhoIsoWithEACut_C1 + 0.0038 # absPFPhoIsoWithEACut_C2 + ) + +# Tight working point Barrel and Endcap +idName = "cutBasedPhotonID-Fall17-94X-V1-tight" +WP_Tight_EB = WorkingPoint_V2( + idName , # idName + 0.022 , # hOverECut + 0.0099 , # full5x5_SigmaIEtaIEtaCut +# Isolation cuts are generally absIso < C1 + pt*C2, except for NeuHad is < C1 + pt*C2 + pt*pt*C3 + 0.101 , # absPFChaHadIsoWithEACut_C1 + 0.00 , # absPFChaHadIsoWithEACut_C2 + 0.137 , # absPFNeuHadIsoWithEACut_C1 + 0.0127 , # absPFNeuHadIsoWithEACut_C2 + 0.000026 , # absPFNeuHadIsowithEACut_C3 + 2.308 , # absPFPhoIsoWithEACut_C1 + 0.0038 # absPFPhoIsoWithEACut_C2 + ) + +WP_Tight_EE = WorkingPoint_V2( + idName , #idName + 0.021 , # hOverECut + 0.0267 , # full5x5_SigmaIEtaIEtaCut +# Isolation cuts are generally absIso < C1 + pt*C2, except for NeuHad is < C1 + pt*C2 + pt*pt*C3 + 0.134 , # absPFChaHadIsoWithEACut_C1 + 0.00 , # absPFChaHadIsoWithEACut_C2 + 1.615 , # absPFNeuHadIsoWithEACut_C1 + 0.0117 , # absPFNeuHadIsoWithEACut_C2 + 0.000025 , # absPFNeuHadIsowithEACut_C3 + 3.107 , # absPFPhoIsoWithEACut_C1 + 0.0038 # absPFPhoIsoWithEACut_C2 + ) + + +# Second, define where to find the precomputed isolations and what effective +# areas to use for pile-up correction +isoInputs = IsolationCutInputs( + # chHadIsolationMapName + 'photonIDValueMapProducer:phoChargedIsolation' , + # chHadIsolationEffAreas + "RecoEgamma/PhotonIdentification/data/Fall17/effAreaPhotons_cone03_pfChargedHadrons_90percentBased.txt", + # neuHadIsolationMapName + 'photonIDValueMapProducer:phoNeutralHadronIsolation' , + # neuHadIsolationEffAreas + "RecoEgamma/PhotonIdentification/data/Fall17/effAreaPhotons_cone03_pfNeutralHadrons_90percentBased.txt" , + # phoIsolationMapName + "photonIDValueMapProducer:phoPhotonIsolation" , + # phoIsolationEffAreas + "RecoEgamma/PhotonIdentification/data/Fall17/effAreaPhotons_cone03_pfPhotons_90percentBased.txt" +) + +# +# Finally, set up VID configuration for all cuts +# +cutBasedPhotonID_Fall17_94X_V1_loose = configureVIDCutBasedPhoID_V5 ( WP_Loose_EB, WP_Loose_EE, isoInputs) +cutBasedPhotonID_Fall17_94X_V1_medium = configureVIDCutBasedPhoID_V5 ( WP_Medium_EB, WP_Medium_EE, isoInputs) +cutBasedPhotonID_Fall17_94X_V1_tight = configureVIDCutBasedPhoID_V5 ( WP_Tight_EB, WP_Tight_EE, isoInputs) + +## 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(cutBasedPhotonID_Fall17_94X_V1_loose.idName, + '08547098f52eb608b545953f02583c3f') +central_id_registry.register(cutBasedPhotonID_Fall17_94X_V1_medium.idName, + 'fb58ccd713d6be1f86f1d2e48c69e401') +central_id_registry.register(cutBasedPhotonID_Fall17_94X_V1_tight.idName, + '296da1cdbf6f35a99287c5a527472ed3') +cutBasedPhotonID_Fall17_94X_V1_loose.isPOGApproved = cms.untracked.bool(True) +cutBasedPhotonID_Fall17_94X_V1_medium.isPOGApproved = cms.untracked.bool(True) +cutBasedPhotonID_Fall17_94X_V1_tight.isPOGApproved = cms.untracked.bool(True) diff --git a/RecoEgamma/PhotonIdentification/python/Identification/mvaPhotonID_Fall17_94X_V1_cff.py b/RecoEgamma/PhotonIdentification/python/Identification/mvaPhotonID_Fall17_94X_V1_cff.py new file mode 100644 index 0000000000000..64b9d0a0c75d2 --- /dev/null +++ b/RecoEgamma/PhotonIdentification/python/Identification/mvaPhotonID_Fall17_94X_V1_cff.py @@ -0,0 +1,126 @@ +from PhysicsTools.SelectorUtils.centralIDRegistry import central_id_registry + +import FWCore.ParameterSet.Config as cms + +# +# In this file we define the locations of the MVA weights, cuts on the MVA values +# for specific working points, and configure those cuts in VID +# + +# +# The following MVA is derived for Fall17 samples for non-triggering photons. +# See more documentation in these presentations: +# https://indico.cern.ch/event/662751/contributions/2778043/attachments/1562017/2459674/EGamma_WorkShop_21.11.17_Debabrata.pdf +# + +# This MVA implementation class name +mvaFall17v1ClassName = "PhotonMVAEstimatorRunIIFall17" +# The tag is an extra string attached to the names of the products +# such as ValueMaps that needs to distinguish cases when the same MVA estimator +# class is used with different tuning/weights +mvaTag = "v1" + +# There are 2 categories in this MVA. They have to be configured in this strict order +# (cuts and weight files order): +# 0 barrel photons +# 1 endcap photons + +mvaRunIIFall17WeightFiles_V1 = cms.vstring( + "RecoEgamma/PhotonIdentification/data/Fall17/HggPhoId_92X_barrel_BDT.weights.xml", + "RecoEgamma/PhotonIdentification/data/Fall17/HggPhoId_92X_endcap_BDT.weights.xml" + ) + +effAreasPath_pho = "RecoEgamma/PhotonIdentification/data/Fall17/effAreaPhotons_cone03_pfPhotons_90percentBased_TrueVtx.txt" + +# Load some common definitions for MVA machinery +from RecoEgamma.PhotonIdentification.Identification.mvaPhotonID_tools \ + import ( PhoMVA_2Categories_WP, + configureVIDMVAPhoID_V1 ) + +# The locatoins of value maps with the actual MVA values and categories +# for all particles. +# The names for the maps are ":Values" +# and ":Categories" +mvaProducerModuleLabel = "photonMVAValueMapProducer" +mvaValueMapName = mvaProducerModuleLabel + ":" + mvaFall17v1ClassName + mvaTag + "Values" +mvaCategoriesMapName = mvaProducerModuleLabel + ":" + mvaFall17v1ClassName + mvaTag + "Categories" + +# The working point for this MVA that is expected to have about 90% signal +# efficiency in each category for photons with pt>30 GeV (somewhat lower +# for lower pt photons). +idName = "mvaPhoID-RunIIFall17-v1-wp90" +MVA_WP90 = PhoMVA_2Categories_WP( + idName = idName, + mvaValueMapName = mvaValueMapName, # map with MVA values for all particles + mvaCategoriesMapName = mvaCategoriesMapName, # map with category index for all particles + cutCategory0 = 0.27, # EB new val : sig eff = 90% , bkg eff = ? + cutCategory1 = 0.14 # EE new val : sig eff = 90% , bkg eff = ? + ) + +# The working point for this MVA that is expected to have about 90% signal +# efficiency in each category for photons with pt>30 GeV (somewhat lower +# for lower pt photons). +idName = "mvaPhoID-RunIIFall17-v1-wp80" +MVA_WP80 = PhoMVA_2Categories_WP( + idName = idName, + mvaValueMapName = mvaValueMapName, # map with MVA values for all particles + mvaCategoriesMapName = mvaCategoriesMapName, # map with category index for all particles + cutCategory0 = 0.67, # EB new val : sig eff = 80% , bkg eff = ? + cutCategory1 = 0.54 # EE new val : sig eff = 80% , bkg eff = ? + ) + +# +# Finally, set up VID configuration for all cuts +# + +# Create the PSet that will be fed to the MVA value map producer +mvaPhoID_RunIIFall17_v1_producer_config = cms.PSet( + mvaName = cms.string(mvaFall17v1ClassName), + mvaTag = cms.string(mvaTag), + weightFileNames = mvaRunIIFall17WeightFiles_V1, + # + # All the event content needed for this MVA implementation follows + # + # All the value maps: these are expected to be produced by the + # PhotonIDValueMapProducer running upstream + # + phoChargedIsolation = cms.InputTag("photonIDValueMapProducer:phoChargedIsolation"), + phoPhotonIsolation = cms.InputTag("photonIDValueMapProducer:phoPhotonIsolation"), + phoWorstChargedIsolation = cms.InputTag("photonIDValueMapProducer:phoWorstChargedIsolationWithPVConstraint"), + # + # Original event content: pileup in this case + # + rho = cms.InputTag("fixedGridRhoAll"), # As used by Hgg and by developer of this ID + # In this MVA for endcap the corrected photon isolation is defined as + # iso = max( photon_isolation_raw - rho*effArea - coeff*pt, cutoff) + # as discussed in the indico presentations listed in the beginning of this file. + # + effAreasConfigFile = cms.FileInPath(effAreasPath_pho), + # The coefficients "coeff" for the formula above for linear pt scaling correction + # the first value is for EB, the second is for EE + # NOTE: even though the EB coefficient is provided, it is not presently used in the MVA. + # For EB, the uncorrected raw photon isolation is used instead. + phoIsoPtScalingCoeff = cms.vdouble(0.0035,0.0040) + # The cutoff for the formula above + # phoIsoCutoff = cms.double(2.5) + ) + +# Create the VPset's for VID cuts +mvaPhoID_RunIIFall17_v1_wp90 = configureVIDMVAPhoID_V1( MVA_WP90 ) +mvaPhoID_RunIIFall17_v1_wp80 = configureVIDMVAPhoID_V1( MVA_WP80 ) + +# 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 "calculateIdMD5 +# 3) update the MD5 sum strings below and uncomment the lines again. +# + + +central_id_registry.register( mvaPhoID_RunIIFall17_v1_wp90.idName, + '834dd792692b6a62786bd1caa6b53a68') +central_id_registry.register( mvaPhoID_RunIIFall17_v1_wp80.idName, + '7e48c47329d7d1eb889100ed03a02ba9') + +mvaPhoID_RunIIFall17_v1_wp90.isPOGApproved = cms.untracked.bool(True) +mvaPhoID_RunIIFall17_v1_wp80.isPOGApproved = cms.untracked.bool(True) diff --git a/RecoEgamma/PhotonIdentification/python/Identification/mvaPhotonID_Fall17_94X_V1p1_cff.py b/RecoEgamma/PhotonIdentification/python/Identification/mvaPhotonID_Fall17_94X_V1p1_cff.py new file mode 100644 index 0000000000000..2a992f446dd2e --- /dev/null +++ b/RecoEgamma/PhotonIdentification/python/Identification/mvaPhotonID_Fall17_94X_V1p1_cff.py @@ -0,0 +1,126 @@ +from PhysicsTools.SelectorUtils.centralIDRegistry import central_id_registry + +import FWCore.ParameterSet.Config as cms + +# +# In this file we define the locations of the MVA weights, cuts on the MVA values +# for specific working points, and configure those cuts in VID +# + +# +# The following MVA is derived for Fall17 samples for non-triggering photons. +# See more documentation in these presentations: +# https://indico.cern.ch/event/662751/contributions/2778043/attachments/1562017/2459674/EGamma_WorkShop_21.11.17_Debabrata.pdf +# + +# This MVA implementation class name +mvaFall17v1p1ClassName = "PhotonMVAEstimatorRunIIFall17" +# The tag is an extra string attached to the names of the products +# such as ValueMaps that needs to distinguish cases when the same MVA estimator +# class is used with different tuning/weights +mvaTag = "v1p1" + +# There are 2 categories in this MVA. They have to be configured in this strict order +# (cuts and weight files order): +# 0 barrel photons +# 1 endcap photons + +mvaRunIIFall17WeightFiles_V1p1 = cms.vstring( + "RecoEgamma/PhotonIdentification/data/Fall17/HggPhoId_92X_barrel_BDT.weights.xml", + "RecoEgamma/PhotonIdentification/data/Fall17/HggPhoId_92X_endcap_BDT.weights.xml" + ) + +effAreasPath_pho = "RecoEgamma/PhotonIdentification/data/Fall17/effAreaPhotons_cone03_pfPhotons_90percentBased_TrueVtx.txt" + +# Load some common definitions for MVA machinery +from RecoEgamma.PhotonIdentification.Identification.mvaPhotonID_tools \ + import ( PhoMVA_2Categories_WP, + configureVIDMVAPhoID_V1 ) + +# The locatoins of value maps with the actual MVA values and categories +# for all particles. +# The names for the maps are ":Values" +# and ":Categories" +mvaProducerModuleLabel = "photonMVAValueMapProducer" +mvaValueMapName = mvaProducerModuleLabel + ":" + mvaFall17v1p1ClassName + mvaTag + "Values" +mvaCategoriesMapName = mvaProducerModuleLabel + ":" + mvaFall17v1p1ClassName + mvaTag + "Categories" + +# The working point for this MVA that is expected to have about 90% signal +# efficiency in each category for photons with pt>30 GeV (somewhat lower +# for lower pt photons). +idName = "mvaPhoID-RunIIFall17-v1p1-wp90" +MVA_WP90 = PhoMVA_2Categories_WP( + idName = idName, + mvaValueMapName = mvaValueMapName, # map with MVA values for all particles + mvaCategoriesMapName = mvaCategoriesMapName, # map with category index for all particles + cutCategory0 = 0.27, # EB new val : sig eff = 90% , bkg eff = ? + cutCategory1 = 0.14 # EE new val : sig eff = 90% , bkg eff = ? + ) + +# The working point for this MVA that is expected to have about 90% signal +# efficiency in each category for photons with pt>30 GeV (somewhat lower +# for lower pt photons). +idName = "mvaPhoID-RunIIFall17-v1p1-wp80" +MVA_WP80 = PhoMVA_2Categories_WP( + idName = idName, + mvaValueMapName = mvaValueMapName, # map with MVA values for all particles + mvaCategoriesMapName = mvaCategoriesMapName, # map with category index for all particles + cutCategory0 = 0.67, # EB new val : sig eff = 80% , bkg eff = ? + cutCategory1 = 0.54 # EE new val : sig eff = 80% , bkg eff = ? + ) + +# +# Finally, set up VID configuration for all cuts +# + +# Create the PSet that will be fed to the MVA value map producer +mvaPhoID_RunIIFall17_v1p1_producer_config = cms.PSet( + mvaName = cms.string(mvaFall17v1p1ClassName), + mvaTag = cms.string(mvaTag), + weightFileNames = mvaRunIIFall17WeightFiles_V1p1, + # + # All the event content needed for this MVA implementation follows + # + # All the value maps: these are expected to be produced by the + # PhotonIDValueMapProducer running upstream + # + phoChargedIsolation = cms.InputTag("photonIDValueMapProducer:phoChargedIsolation"), + phoPhotonIsolation = cms.InputTag("photonIDValueMapProducer:phoPhotonIsolation"), + phoWorstChargedIsolation = cms.InputTag("photonIDValueMapProducer:phoWorstChargedIsolation"), + # + # Original event content: pileup in this case + # + rho = cms.InputTag("fixedGridRhoAll"), # As used by Hgg and by developer of this ID + # In this MVA for endcap the corrected photon isolation is defined as + # iso = max( photon_isolation_raw - rho*effArea - coeff*pt, cutoff) + # as discussed in the indico presentations listed in the beginning of this file. + # + effAreasConfigFile = cms.FileInPath(effAreasPath_pho), + # The coefficients "coeff" for the formula above for linear pt scaling correction + # the first value is for EB, the second is for EE + # NOTE: even though the EB coefficient is provided, it is not presently used in the MVA. + # For EB, the uncorrected raw photon isolation is used instead. + phoIsoPtScalingCoeff = cms.vdouble(0.0035,0.0040) + # The cutoff for the formula above + # phoIsoCutoff = cms.double(2.5) + ) + +# Create the VPset's for VID cuts +mvaPhoID_RunIIFall17_v1p1_wp90 = configureVIDMVAPhoID_V1( MVA_WP90 ) +mvaPhoID_RunIIFall17_v1p1_wp80 = configureVIDMVAPhoID_V1( MVA_WP80 ) + +# 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 "calculateIdMD5 +# 3) update the MD5 sum strings below and uncomment the lines again. +# + + +central_id_registry.register( mvaPhoID_RunIIFall17_v1p1_wp90.idName, + '1120f91d15f68bf61b5f08958bf4f435') +central_id_registry.register( mvaPhoID_RunIIFall17_v1p1_wp80.idName, + '56138c4a3ac3c0bffc7f01c187063102') + +mvaPhoID_RunIIFall17_v1p1_wp90.isPOGApproved = cms.untracked.bool(True) +mvaPhoID_RunIIFall17_v1p1_wp80.isPOGApproved = cms.untracked.bool(True) diff --git a/RecoEgamma/PhotonIdentification/python/PhotonMVAValueMapProducer_cfi.py b/RecoEgamma/PhotonIdentification/python/PhotonMVAValueMapProducer_cfi.py index 1420ed2631767..c22de64d0a3f4 100644 --- a/RecoEgamma/PhotonIdentification/python/PhotonMVAValueMapProducer_cfi.py +++ b/RecoEgamma/PhotonIdentification/python/PhotonMVAValueMapProducer_cfi.py @@ -14,6 +14,12 @@ import mvaPhoID_Spring16_nonTrig_V1_producer_config mvaConfigsForPhoProducer.append( mvaPhoID_Spring16_nonTrig_V1_producer_config ) +from RecoEgamma.PhotonIdentification.Identification.mvaPhotonID_Fall17_94X_V1_cff import * +mvaConfigsForPhoProducer.append( mvaPhoID_RunIIFall17_v1_producer_config ) + +from RecoEgamma.PhotonIdentification.Identification.mvaPhotonID_Fall17_94X_V1p1_cff import * +mvaConfigsForPhoProducer.append( mvaPhoID_RunIIFall17_v1p1_producer_config ) + photonMVAValueMapProducer = cms.EDProducer('PhotonMVAValueMapProducer', # The module automatically detects AOD vs miniAOD, so we configure both # diff --git a/RecoEgamma/PhotonIdentification/test/runtests.sh b/RecoEgamma/PhotonIdentification/test/runtests.sh index 79fdce43ccc65..4028801c4ad13 100755 --- a/RecoEgamma/PhotonIdentification/test/runtests.sh +++ b/RecoEgamma/PhotonIdentification/test/runtests.sh @@ -1,10 +1,14 @@ function die { echo $1: status $2 ; exit $2; } ids_to_test=( -"RecoEgamma.PhotonIdentification.Identification.cutBasedPhotonID_Spring15_25ns_V1_cff" -"RecoEgamma.PhotonIdentification.Identification.cutBasedPhotonID_Spring15_50ns_V1_cff" -"RecoEgamma.PhotonIdentification.Identification.mvaPhotonID_Spring15_25ns_nonTrig_V2p1_cff" -"RecoEgamma.PhotonIdentification.Identification.mvaPhotonID_Spring15_50ns_nonTrig_V2p1_cff" + "RecoEgamma.PhotonIdentification.Identification.cutBasedPhotonID_Spring15_25ns_V1_cff" + "RecoEgamma.PhotonIdentification.Identification.cutBasedPhotonID_Spring15_50ns_V1_cff" + "RecoEgamma.PhotonIdentification.Identification.mvaPhotonID_Spring15_25ns_nonTrig_V2p1_cff" + "RecoEgamma.PhotonIdentification.Identification.mvaPhotonID_Spring15_50ns_nonTrig_V2p1_cff" + 'RecoEgamma.PhotonIdentification.Identification.cutBasedPhotonID_Spring16_V2p2_cff' + 'RecoEgamma.PhotonIdentification.Identification.mvaPhotonID_Spring16_nonTrig_V1_cff' + 'RecoEgamma.PhotonIdentification.Identification.cutBasedPhotonID_Fall17_94X_V1_TrueVtx_cff' + 'RecoEgamma.PhotonIdentification.Identification.mvaPhotonID_Fall17_94X_V1_cff' ) for id_set in "${ids_to_test[@]}"; do