diff --git a/DataFormats/EgammaCandidates/interface/GsfElectron.h b/DataFormats/EgammaCandidates/interface/GsfElectron.h index b81bdf23d1974..22b6c65bd79ef 100644 --- a/DataFormats/EgammaCandidates/interface/GsfElectron.h +++ b/DataFormats/EgammaCandidates/interface/GsfElectron.h @@ -28,8 +28,8 @@ namespace reco { * Renamed from PixelMatchGsfElectron. * Originally adapted from the TRecElectron class in ORCA. * - * \author Claude Charlot - Laboratoire Leprince-Ringuet - École polytechnique, CNRS/IN2P3 - * \author David Chamont - Laboratoire Leprince-Ringuet - École polytechnique, CNRS/IN2P3 + * \author Claude Charlot - Laboratoire Leprince-Ringuet - École polytechnique, CNRS/IN2P3 + * \author David Chamont - Laboratoire Leprince-Ringuet - École polytechnique, CNRS/IN2P3 * ****************************************************************************/ @@ -701,12 +701,28 @@ namespace reco { etOutsideMustache(-std::numeric_limits::max()) {} }; + static constexpr float mvaPlaceholder = -999999999.; + struct MvaOutput { int status; // see PFCandidateElectronExtra::StatusFlag float mva_Isolated; float mva_e_pi; float mvaByPassForIsolated; // complementary MVA used in preselection - MvaOutput() : status(-1), mva_Isolated(-999999999.), mva_e_pi(-999999999.), mvaByPassForIsolated(-999999999.) {} + float dnn_e_sigIsolated; + float dnn_e_sigNonIsolated; + float dnn_e_bkgNonIsolated; + float dnn_e_bkgTau; + float dnn_e_bkgPhoton; + MvaOutput() + : status(-1), + mva_Isolated(mvaPlaceholder), + mva_e_pi(mvaPlaceholder), + mvaByPassForIsolated(mvaPlaceholder), + dnn_e_sigIsolated(mvaPlaceholder), + dnn_e_sigNonIsolated(mvaPlaceholder), + dnn_e_bkgNonIsolated(mvaPlaceholder), + dnn_e_bkgTau(mvaPlaceholder), + dnn_e_bkgPhoton(mvaPlaceholder) {} }; // accessors @@ -726,6 +742,11 @@ namespace reco { // for backward compatibility float mva_Isolated() const { return mvaOutput_.mva_Isolated; } float mva_e_pi() const { return mvaOutput_.mva_e_pi; } + float dnn_signal_Isolated() const { return mvaOutput_.dnn_e_sigIsolated; } + float dnn_signal_nonIsolated() const { return mvaOutput_.dnn_e_sigNonIsolated; } + float dnn_bkg_nonIsolated() const { return mvaOutput_.dnn_e_bkgNonIsolated; } + float dnn_bkg_Tau() const { return mvaOutput_.dnn_e_bkgTau; } + float dnn_bkg_Photon() const { return mvaOutput_.dnn_e_bkgPhoton; } private: PflowIsolationVariables pfIso_; diff --git a/DataFormats/EgammaCandidates/interface/Photon.h b/DataFormats/EgammaCandidates/interface/Photon.h index 071f0329a03fb..f1daf95615eb4 100644 --- a/DataFormats/EgammaCandidates/interface/Photon.h +++ b/DataFormats/EgammaCandidates/interface/Photon.h @@ -556,25 +556,23 @@ namespace reco { /// Set Particle Flow Isolation variables void setPflowIsolationVariables(const PflowIsolationVariables& pfisol) { pfIsolation_ = pfisol; } + static constexpr float mvaPlaceholder = -999999999.; + struct PflowIDVariables { int nClusterOutsideMustache; float etOutsideMustache; float mva; + float dnn; PflowIDVariables() - : - - nClusterOutsideMustache(-1), - etOutsideMustache(-999999999.), - mva(-999999999.) - - {} + : nClusterOutsideMustache(-1), etOutsideMustache(mvaPlaceholder), mva(mvaPlaceholder), dnn(mvaPlaceholder) {} }; // getters int nClusterOutsideMustache() const { return pfID_.nClusterOutsideMustache; } float etOutsideMustache() const { return pfID_.etOutsideMustache; } float pfMVA() const { return pfID_.mva; } + float pfDNN() const { return pfID_.dnn; } // setters void setPflowIDVariables(const PflowIDVariables& pfid) { pfID_ = pfid; } diff --git a/DataFormats/EgammaCandidates/src/classes_def.xml b/DataFormats/EgammaCandidates/src/classes_def.xml index 7b93b9863aa35..9432329160786 100644 --- a/DataFormats/EgammaCandidates/src/classes_def.xml +++ b/DataFormats/EgammaCandidates/src/classes_def.xml @@ -69,7 +69,8 @@ - + + @@ -221,7 +222,8 @@ - + + diff --git a/DataFormats/ParticleFlowCandidate/interface/PFCandidate.h b/DataFormats/ParticleFlowCandidate/interface/PFCandidate.h index 31d43fa4c9087..9ca9bce8832c5 100644 --- a/DataFormats/ParticleFlowCandidate/interface/PFCandidate.h +++ b/DataFormats/ParticleFlowCandidate/interface/PFCandidate.h @@ -343,6 +343,31 @@ namespace reco { /// set mva for neutral hadron - gamma discrimination void set_mva_gamma_nh(float mva) { mva_gamma_nh_ = mva; } + // set DNN for electron PFID + // mva for ele PFID DNN sigIsolated class + float dnn_e_sigIsolated() const { return dnn_e_sigIsolated_; } + void set_dnn_e_sigIsolated(float mva) { dnn_e_sigIsolated_ = mva; } + + // mva for ele PFID DNN sigNonIsolated class + float dnn_e_sigNonIsolated() const { return dnn_e_sigNonIsolated_; } + void set_dnn_e_sigNonIsolated(float mva) { dnn_e_sigNonIsolated_ = mva; } + + // mva for ele PFID DNN bkgNonIsolated class + float dnn_e_bkgNonIsolated() const { return dnn_e_bkgNonIsolated_; } + void set_dnn_e_bkgNonIsolated(float mva) { dnn_e_bkgNonIsolated_ = mva; } + + // mva for ele PFID DNN bkgTau class + float dnn_e_bkgTau() const { return dnn_e_bkgTau_; } + void set_dnn_e_bkgTau(float mva) { dnn_e_bkgTau_ = mva; } + + // mva for ele PFID DNN bkgPhoton class + float dnn_e_bkgPhoton() const { return dnn_e_bkgPhoton_; } + void set_dnn_e_bkgPhoton(float mva) { dnn_e_bkgPhoton_ = mva; } + + // set DNN for gamma PFID + float dnn_gamma() const { return dnn_gamma_; } + void set_dnn_gamma(float mva) { dnn_gamma_ = mva; } + /// mva for neutral hadron - gamma discrimination float mva_gamma_nh() const { return mva_gamma_nh_; } @@ -383,7 +408,7 @@ namespace reco { const ElementsInBlocks& elementsInBlocks() const; - static const float bigMva_; + static constexpr float bigMva_ = -999.; friend std::ostream& operator<<(std::ostream& out, const PFCandidate& c); @@ -492,6 +517,24 @@ namespace reco { /// mva for neutral hadron - gamma discrimination float mva_gamma_nh_; + /// DNN for electron PFid: isolated signal + float dnn_e_sigIsolated_; + + /// DNN for electron PFid: non-isolated signal + float dnn_e_sigNonIsolated_; + + /// DNN for electron PFid: non-isolated bkg + float dnn_e_bkgNonIsolated_; + + /// DNN for electron PFid: tau bkg + float dnn_e_bkgTau_; + + /// DNN for electron PFid: photon bkg + float dnn_e_bkgPhoton_; + + // DNN for gamma PFid + float dnn_gamma_; + /// position at ECAL entrance, from the PFRecTrack math::XYZPointF positionAtECALEntrance_; diff --git a/DataFormats/ParticleFlowCandidate/src/PFCandidate.cc b/DataFormats/ParticleFlowCandidate/src/PFCandidate.cc index c178c042e604b..632579a364f6d 100644 --- a/DataFormats/ParticleFlowCandidate/src/PFCandidate.cc +++ b/DataFormats/ParticleFlowCandidate/src/PFCandidate.cc @@ -18,8 +18,6 @@ using namespace reco; using namespace std; -const float PFCandidate::bigMva_ = -999.; - #include "DataFormats/ParticleFlowCandidate/src/CountBits.h" PFCandidate::PFCandidate() @@ -35,13 +33,18 @@ PFCandidate::PFCandidate() flags_(0), deltaP_(0.), vertexType_(kCandVertex), - mva_Isolated_(bigMva_), - mva_e_pi_(bigMva_), - mva_e_mu_(bigMva_), - mva_pi_mu_(bigMva_), - mva_nothing_gamma_(bigMva_), - mva_nothing_nh_(bigMva_), - mva_gamma_nh_(bigMva_), + mva_Isolated_(PFCandidate::bigMva_), + mva_e_pi_(PFCandidate::bigMva_), + mva_e_mu_(PFCandidate::bigMva_), + mva_pi_mu_(PFCandidate::bigMva_), + mva_nothing_gamma_(PFCandidate::bigMva_), + mva_nothing_nh_(PFCandidate::bigMva_), + mva_gamma_nh_(PFCandidate::bigMva_), + dnn_e_sigIsolated_(PFCandidate::bigMva_), + dnn_e_sigNonIsolated_(PFCandidate::bigMva_), + dnn_e_bkgNonIsolated_(PFCandidate::bigMva_), + dnn_e_bkgTau_(PFCandidate::bigMva_), + dnn_e_bkgPhoton_(PFCandidate::bigMva_), getter_(nullptr), storedRefsBitPattern_(0), time_(0.f), @@ -74,13 +77,18 @@ PFCandidate::PFCandidate(Charge charge, const LorentzVector& p4, ParticleType pa flags_(0), deltaP_(0.), vertexType_(kCandVertex), - mva_Isolated_(bigMva_), - mva_e_pi_(bigMva_), - mva_e_mu_(bigMva_), - mva_pi_mu_(bigMva_), - mva_nothing_gamma_(bigMva_), - mva_nothing_nh_(bigMva_), - mva_gamma_nh_(bigMva_), + mva_Isolated_(PFCandidate::bigMva_), + mva_e_pi_(PFCandidate::bigMva_), + mva_e_mu_(PFCandidate::bigMva_), + mva_pi_mu_(PFCandidate::bigMva_), + mva_nothing_gamma_(PFCandidate::bigMva_), + mva_nothing_nh_(PFCandidate::bigMva_), + mva_gamma_nh_(PFCandidate::bigMva_), + dnn_e_sigIsolated_(PFCandidate::bigMva_), + dnn_e_sigNonIsolated_(PFCandidate::bigMva_), + dnn_e_bkgNonIsolated_(PFCandidate::bigMva_), + dnn_e_bkgTau_(PFCandidate::bigMva_), + dnn_e_bkgPhoton_(PFCandidate::bigMva_), getter_(nullptr), storedRefsBitPattern_(0), time_(0.f), @@ -137,6 +145,11 @@ PFCandidate::PFCandidate(PFCandidate const& iOther) mva_nothing_gamma_(iOther.mva_nothing_gamma_), mva_nothing_nh_(iOther.mva_nothing_nh_), mva_gamma_nh_(iOther.mva_gamma_nh_), + dnn_e_sigIsolated_(iOther.dnn_e_sigIsolated_), + dnn_e_sigNonIsolated_(iOther.dnn_e_sigNonIsolated_), + dnn_e_bkgNonIsolated_(iOther.dnn_e_bkgNonIsolated_), + dnn_e_bkgTau_(iOther.dnn_e_bkgTau_), + dnn_e_bkgPhoton_(iOther.dnn_e_bkgPhoton_), positionAtECALEntrance_(iOther.positionAtECALEntrance_), getter_(iOther.getter_), storedRefsBitPattern_(iOther.storedRefsBitPattern_), @@ -181,6 +194,11 @@ PFCandidate& PFCandidate::operator=(PFCandidate const& iOther) { mva_nothing_gamma_ = iOther.mva_nothing_gamma_; mva_nothing_nh_ = iOther.mva_nothing_nh_; mva_gamma_nh_ = iOther.mva_gamma_nh_; + dnn_e_sigIsolated_ = iOther.dnn_e_sigIsolated_; + dnn_e_sigNonIsolated_ = iOther.dnn_e_sigNonIsolated_; + dnn_e_bkgNonIsolated_ = iOther.dnn_e_bkgNonIsolated_; + dnn_e_bkgTau_ = iOther.dnn_e_bkgTau_; + dnn_e_bkgPhoton_ = iOther.dnn_e_bkgPhoton_; positionAtECALEntrance_ = iOther.positionAtECALEntrance_; getter_ = iOther.getter_; storedRefsBitPattern_ = iOther.storedRefsBitPattern_; diff --git a/DataFormats/ParticleFlowCandidate/src/classes_def.xml b/DataFormats/ParticleFlowCandidate/src/classes_def.xml index eca906716f5bc..17fa11e0ed1ee 100644 --- a/DataFormats/ParticleFlowCandidate/src/classes_def.xml +++ b/DataFormats/ParticleFlowCandidate/src/classes_def.xml @@ -1,6 +1,7 @@ - + + @@ -52,7 +53,8 @@ - + + @@ -67,7 +69,8 @@ - + + diff --git a/DataFormats/PatCandidates/src/classes_def_objects.xml b/DataFormats/PatCandidates/src/classes_def_objects.xml index e4a78ebf00e43..2b98aef02aab2 100644 --- a/DataFormats/PatCandidates/src/classes_def_objects.xml +++ b/DataFormats/PatCandidates/src/classes_def_objects.xml @@ -314,7 +314,8 @@ - + + diff --git a/RecoEgamma/EgammaElectronAlgos/interface/GsfElectronAlgo.h b/RecoEgamma/EgammaElectronAlgos/interface/GsfElectronAlgo.h index 585edada3e9db..65c27221a9133 100644 --- a/RecoEgamma/EgammaElectronAlgos/interface/GsfElectronAlgo.h +++ b/RecoEgamma/EgammaElectronAlgos/interface/GsfElectronAlgo.h @@ -37,6 +37,7 @@ #include "RecoEgamma/EgammaIsolationAlgos/interface/EleTkIsolFromCands.h" #include "RecoEgamma/ElectronIdentification/interface/ElectronMVAEstimator.h" #include "RecoEgamma/ElectronIdentification/interface/SoftElectronMVAEstimator.h" +#include "RecoEgamma/ElectronIdentification/interface/ElectronDNNEstimator.h" #include "RecoLocalCalo/EcalRecAlgos/interface/EcalSeverityLevelAlgo.h" #include "RecoLocalCalo/EcalRecAlgos/interface/EcalSeverityLevelAlgoRcd.h" #include "TrackingTools/GsfTools/interface/MultiTrajectoryStateMode.h" @@ -47,6 +48,8 @@ #include "RecoEgamma/EgammaElectronAlgos/interface/ConversionFinder.h" #include "CondFormats/EcalObjects/interface/EcalPFRecHitThresholds.h" #include "CondFormats/DataRecord/interface/EcalPFRecHitThresholdsRcd.h" +#include "RecoEgamma/EgammaIsolationAlgos/interface/EcalPFClusterIsolation.h" +#include "RecoEgamma/EgammaIsolationAlgos/interface/HcalPFClusterIsolation.h" class GsfElectronAlgo { public: @@ -55,6 +58,7 @@ class GsfElectronAlgo { HeavyObjectCache(const edm::ParameterSet&); std::unique_ptr sElectronMVAEstimator; std::unique_ptr iElectronMVAEstimator; + std::unique_ptr iElectronDNNEstimator; }; struct Tokens { @@ -69,6 +73,11 @@ class GsfElectronAlgo { edm::EDGetTokenT beamSpotTag; edm::EDGetTokenT vtxCollectionTag; edm::EDGetTokenT conversions; + + edm::EDGetTokenT pfClusterProducer; + edm::EDGetTokenT pfClusterProducerHCAL; + edm::EDGetTokenT pfClusterProducerHFEM; + edm::EDGetTokenT pfClusterProducerHFHAD; }; struct StrategyConfiguration { @@ -95,6 +104,8 @@ class GsfElectronAlgo { //heavy ion in 2015 has no conversions and so cant fill conv vtx fit prob so this bool //stops it from being filled bool fillConvVtxFitProb; + // Compute PFcluster isolation for egamma PFID DNN + bool computePfClusterIso; }; struct CutsConfiguration { @@ -175,12 +186,33 @@ class GsfElectronAlgo { bool useNumCrystals; }; + struct PFClusterIsolationConfiguration { + double ecaldrMax; + double ecaldrVetoBarrel; + double ecaldrVetoEndcap; + double ecaletaStripBarrel; + double ecaletaStripEndcap; + double ecalenergyBarrel; + double ecalenergyEndcap; + + bool useHF; + double hcaldrMax; + double hcaldrVetoBarrel; + double hcaldrVetoEndcap; + double hcaletaStripBarrel; + double hcaletaStripEndcap; + double hcalenergyBarrel; + double hcalenergyEndcap; + bool hcaluseEt; + }; + GsfElectronAlgo(const Tokens&, const StrategyConfiguration&, const CutsConfiguration& cutsCfg, const ElectronHcalHelper::Configuration& hcalCone, const ElectronHcalHelper::Configuration& hcalBc, const IsolationConfiguration&, + const PFClusterIsolationConfiguration&, const EcalRecHitsConfiguration&, std::unique_ptr&& crackCorrectionFunction, const RegressionHelper::Configuration& regCfg, @@ -204,6 +236,7 @@ class GsfElectronAlgo { const StrategyConfiguration strategy; const CutsConfiguration cuts; const IsolationConfiguration iso; + const PFClusterIsolationConfiguration pfiso; const EcalRecHitsConfiguration recHits; }; @@ -263,6 +296,12 @@ class GsfElectronAlgo { ElectronHcalHelper hcalHelperBc_; std::unique_ptr crackCorrectionFunction_; RegressionHelper regHelper_; + + // Algos for PfCluster Isolation + typedef EcalPFClusterIsolation ElectronEcalPFClusterIsolation; + std::unique_ptr ecalisoAlgo_; + typedef HcalPFClusterIsolation ElectronHcalPFClusterIsolation; + std::unique_ptr hcalisoAlgo_; }; #endif // GsfElectronAlgo_H diff --git a/RecoEgamma/EgammaElectronAlgos/src/GsfElectronAlgo.cc b/RecoEgamma/EgammaElectronAlgos/src/GsfElectronAlgo.cc index 1562e13a0a2f5..7613eb063fd9f 100644 --- a/RecoEgamma/EgammaElectronAlgos/src/GsfElectronAlgo.cc +++ b/RecoEgamma/EgammaElectronAlgos/src/GsfElectronAlgo.cc @@ -45,6 +45,20 @@ GsfElectronAlgo::HeavyObjectCache::HeavyObjectCache(const edm::ParameterSet& con ElectronMVAEstimator::Configuration iconfig; iconfig.vweightsfiles = conf.getParameter>("ElecMVAFilesString"); iElectronMVAEstimator = std::make_unique(iconfig); + + // Here we will have to load the DNN PFID if present in the config + egammaTools::DNNConfiguration dconfig; + const auto& pset_dnn = conf.getParameter("EleDNNPFid"); + const bool dnnEnabled = pset_dnn.getParameter("enabled"); + if (dnnEnabled) { + dconfig.inputTensorName = pset_dnn.getParameter("inputTensorName"); + dconfig.outputTensorName = pset_dnn.getParameter("outputTensorName"); + dconfig.modelsFiles = pset_dnn.getParameter>("modelsFiles"); + dconfig.scalersFiles = pset_dnn.getParameter>("scalersFiles"); + dconfig.outputDim = pset_dnn.getParameter("outputDim"); + const auto useEBModelInGap = pset_dnn.getParameter("useEBModelInGap"); + iElectronDNNEstimator = std::make_unique(dconfig, useEBModelInGap); + } } //=================================================================== @@ -84,6 +98,10 @@ struct GsfElectronAlgo::EventData { bool originalCtfTrackCollectionRetreived = false; bool originalGsfTrackCollectionRetreived = false; + + //PFCluster Isolation handles + edm::Handle ecalClustersHandle; + std::vector> hcalClustersHandle; }; //=================================================================== @@ -381,6 +399,7 @@ GsfElectronAlgo::GsfElectronAlgo(const Tokens& input, const ElectronHcalHelper::Configuration& hcalCone, const ElectronHcalHelper::Configuration& hcalBc, const IsolationConfiguration& iso, + const PFClusterIsolationConfiguration& pfiso, const EcalRecHitsConfiguration& recHits, std::unique_ptr&& crackCorrectionFunction, const RegressionHelper::Configuration& reg, @@ -389,7 +408,7 @@ GsfElectronAlgo::GsfElectronAlgo(const Tokens& input, const edm::ParameterSet& tkIsolHEEP03, const edm::ParameterSet& tkIsolHEEP04, edm::ConsumesCollector&& cc) - : cfg_{input, strategy, cuts, iso, recHits}, + : cfg_{input, strategy, cuts, iso, pfiso, recHits}, tkIsol03CalcCfg_(tkIsol03), tkIsol04CalcCfg_(tkIsol04), tkIsolHEEP03CalcCfg_(tkIsolHEEP03), @@ -405,7 +424,24 @@ GsfElectronAlgo::GsfElectronAlgo(const Tokens& input, crackCorrectionFunction_{std::forward>(crackCorrectionFunction)}, regHelper_{reg, cfg_.strategy.useEcalRegression, cfg_.strategy.useCombinationRegression, cc} -{} +{ + ///PF ECAL cluster based isolations + ecalisoAlgo_ = std::make_unique(pfiso.ecaldrMax, + pfiso.ecaldrVetoBarrel, + pfiso.ecaldrVetoEndcap, + pfiso.ecaletaStripBarrel, + pfiso.ecaletaStripEndcap, + pfiso.ecalenergyBarrel, + pfiso.ecalenergyEndcap); + hcalisoAlgo_ = std::make_unique(pfiso.hcaldrMax, + pfiso.hcaldrVetoBarrel, + pfiso.hcaldrVetoEndcap, + pfiso.hcaletaStripBarrel, + pfiso.hcaletaStripEndcap, + pfiso.hcalenergyBarrel, + pfiso.hcalenergyEndcap, + pfiso.hcaluseEt); +} void GsfElectronAlgo::checkSetup(const edm::EventSetup& es) { if (cfg_.strategy.useEcalRegression || cfg_.strategy.useCombinationRegression) @@ -434,6 +470,17 @@ GsfElectronAlgo::EventData GsfElectronAlgo::beginEvent(edm::Event const& event, auto barrelRecHits = event.getHandle(cfg_.tokens.barrelRecHitCollection); auto endcapRecHits = event.getHandle(cfg_.tokens.endcapRecHitCollection); + edm::Handle ecalPFClusters; + std::vector> hcalPFClusters; + if (cfg_.strategy.computePfClusterIso) { + ecalPFClusters = event.getHandle(cfg_.tokens.pfClusterProducer); + hcalPFClusters.push_back(event.getHandle(cfg_.tokens.pfClusterProducerHCAL)); + if (cfg_.pfiso.useHF) { + hcalPFClusters.push_back(event.getHandle(cfg_.tokens.pfClusterProducerHFEM)); + hcalPFClusters.push_back(event.getHandle(cfg_.tokens.pfClusterProducerHFHAD)); + } + } + auto ctfTracks = event.getHandle(cfg_.tokens.ctfTracks); EventData eventData{ @@ -582,7 +629,12 @@ GsfElectronAlgo::EventData GsfElectronAlgo::beginEvent(edm::Event const& event, .tkIsolHEEP03Calc = EleTkIsolFromCands(tkIsolHEEP03CalcCfg_, *ctfTracks), .tkIsolHEEP04Calc = EleTkIsolFromCands(tkIsolHEEP04CalcCfg_, *ctfTracks), .originalCtfTracks = {}, - .originalGsfTracks = {}}; + .originalGsfTracks = {}, + + .ecalClustersHandle = ecalPFClusters, + .hcalClustersHandle = hcalPFClusters + + }; eventData.ecalBarrelIsol03.setUseNumCrystals(cfg_.iso.useNumCrystals); eventData.ecalBarrelIsol03.setVetoClustered(cfg_.iso.vetoClustered); @@ -1122,6 +1174,18 @@ void GsfElectronAlgo::createElectron(reco::GsfElectronCollection& electrons, ele.setIsolation03(dr03); ele.setIsolation04(dr04); + //==================================================== + // PFclusters based ISO ! + // Added in CMSSW_12_0_1 at this stage to be able to use them in PF ID DNN + //==================================================== + if (cfg_.strategy.computePfClusterIso) { + reco::GsfElectron::PflowIsolationVariables isoVariables; + isoVariables.sumEcalClusterEt = ecalisoAlgo_->getSum(ele, eventData.ecalClustersHandle); + isoVariables.sumHcalClusterEt = hcalisoAlgo_->getSum(ele, eventData.hcalClustersHandle); + // Other Pfiso variables are initialized at 0 and not used + ele.setPfIsolationVariables(isoVariables); + } + //==================================================== // preselection flag //==================================================== diff --git a/RecoEgamma/EgammaElectronProducers/plugins/BuildFile.xml b/RecoEgamma/EgammaElectronProducers/plugins/BuildFile.xml index 9c8278f27f7f5..8d8cdda990ee7 100644 --- a/RecoEgamma/EgammaElectronProducers/plugins/BuildFile.xml +++ b/RecoEgamma/EgammaElectronProducers/plugins/BuildFile.xml @@ -30,6 +30,7 @@ + diff --git a/RecoEgamma/EgammaElectronProducers/plugins/GsfElectronProducer.cc b/RecoEgamma/EgammaElectronProducers/plugins/GsfElectronProducer.cc index ca978e41a11fd..a1583702d1e05 100644 --- a/RecoEgamma/EgammaElectronProducers/plugins/GsfElectronProducer.cc +++ b/RecoEgamma/EgammaElectronProducers/plugins/GsfElectronProducer.cc @@ -22,6 +22,10 @@ #include "RecoEgamma/EgammaElectronAlgos/interface/ElectronUtilities.h" #include "RecoEgamma/EgammaElectronAlgos/interface/GsfElectronAlgo.h" #include "RecoEcal/EgammaCoreTools/interface/EgammaLocalCovParamDefaults.h" +#include "RecoEgamma/EgammaIsolationAlgos/interface/EcalPFClusterIsolation.h" +#include "RecoEgamma/EgammaIsolationAlgos/interface/HcalPFClusterIsolation.h" + +#include using namespace reco; @@ -29,12 +33,39 @@ namespace { void setMVAOutputs(reco::GsfElectronCollection& electrons, const GsfElectronAlgo::HeavyObjectCache* hoc, - reco::VertexCollection const& vertices) { + reco::VertexCollection const& vertices, + bool dnnPFidEnabled, + const std::vector& tfSessions) { + std::vector mva_outputs(electrons.size()); + size_t iele = 0; for (auto& el : electrons) { GsfElectron::MvaOutput mvaOutput; mvaOutput.mva_e_pi = hoc->sElectronMVAEstimator->mva(el, vertices); mvaOutput.mva_Isolated = hoc->iElectronMVAEstimator->mva(el, vertices.size()); - el.setMvaOutput(mvaOutput); + if (dnnPFidEnabled) { + mva_outputs[iele] = mvaOutput; + } else { + el.setMvaOutput(mvaOutput); + } + iele++; + } + if (dnnPFidEnabled) { + // Here send the list of electrons to the ElectronDNNEstimator and get back the values for all the electrons in one go + LogDebug("GsfElectronProducer") << "Getting DNN PFId for ele"; + const auto& dnn_ele_pfid = hoc->iElectronDNNEstimator->evaluate(electrons, tfSessions); + int jele = 0; + for (auto& el : electrons) { + const auto& values = dnn_ele_pfid[jele]; + // get the previous values + auto& mvaOutput = mva_outputs[jele]; + mvaOutput.dnn_e_sigIsolated = values[0]; + mvaOutput.dnn_e_sigNonIsolated = values[1]; + mvaOutput.dnn_e_bkgNonIsolated = values[2]; + mvaOutput.dnn_e_bkgTau = values[3]; + mvaOutput.dnn_e_bkgPhoton = values[4]; + el.setMvaOutput(mvaOutput); + jele++; + } } } @@ -84,7 +115,9 @@ class GsfElectronProducer : public edm::stream::EDProducer(conf); } - static void globalEndJob(GsfElectronAlgo::HeavyObjectCache const*) {} + void endStream() override; + + static void globalEndJob(GsfElectronAlgo::HeavyObjectCache const*){}; // ------------ method called to produce the data ------------ void produce(edm::Event& event, const edm::EventSetup& setup) override; @@ -116,6 +149,10 @@ class GsfElectronProducer : public edm::stream::EDProducer tfSessions_; }; void GsfElectronProducer::fillDescriptions(edm::ConfigurationDescriptions& descriptions) { @@ -249,6 +286,60 @@ void GsfElectronProducer::fillDescriptions(edm::ConfigurationDescriptions& descr "RecoEgamma/ElectronIdentification/data/TMVA_BDTSoftElectrons_7Feb2014.weights.xml", }); + { + edm::ParameterSetDescription psd1; + psd1.add("enabled", false); + psd1.add("inputTensorName", "FirstLayer_input"); + psd1.add("outputTensorName", "sequential/FinalLayer/Softmax"); + psd1.add("outputDim", 3); // Dimension of output vector + psd1.add>( + "modelsFiles", + {"RecoEgamma/ElectronIdentification/data/Ele_PFID_dnn/lowpT/lowpT_modelDNN.pb", + "RecoEgamma/ElectronIdentification/data/Ele_PFID_dnn/highpTEB/highpTEB_modelDNN.pb", + "RecoEgamma/ElectronIdentification/data/Ele_PFID_dnn/highpTEE/highpTEE_modelDNN.pb"}); + psd1.add>( + "scalersFiles", + {"RecoEgamma/ElectronIdentification/data/Ele_PFID_dnn/lowpT/lowpT_scaler.txt", + "RecoEgamma/ElectronIdentification/data/Ele_PFID_dnn/highpTEB/highpTEB_scaler.txt", + "RecoEgamma/ElectronIdentification/data/Ele_PFID_dnn/highpTEE/highpTEE_scaler.txt"}); + psd1.add("useEBModelInGap", true); + // preselection parameters + desc.add("EleDNNPFid", psd1); + } + + ///For PF cluster isolations + ///ECAL + { + edm::ParameterSetDescription psd0; + psd0.add("pfClusterProducer", edm::InputTag("particleFlowClusterECAL")); + psd0.add("drMax", 0.3); + psd0.add("drVetoBarrel", 0.0); + psd0.add("drVetoEndcap", 0.0); + psd0.add("etaStripBarrel", 0.0); + psd0.add("etaStripEndcap", 0.0); + psd0.add("energyBarrel", 0.0); + psd0.add("energyEndcap", 0.0); + desc.add("pfECALClusIsolCfg", psd0); + } + + ///HCAL + { + edm::ParameterSetDescription psd0; + psd0.add("pfClusterProducerHCAL", edm::InputTag("particleFlowClusterHCAL")); + psd0.add("pfClusterProducerHFEM", edm::InputTag("")); + psd0.add("pfClusterProducerHFHAD", edm::InputTag("")); + psd0.add("useHF", false); + psd0.add("drMax", 0.3); + psd0.add("drVetoBarrel", 0.0); + psd0.add("drVetoEndcap", 0.0); + psd0.add("etaStripBarrel", 0.0); + psd0.add("etaStripEndcap", 0.0); + psd0.add("energyBarrel", 0.0); + psd0.add("energyEndcap", 0.0); + psd0.add("useEt", true); + desc.add("pfHCALClusIsolCfg", psd0); + } + descriptions.add("gsfElectronProducerDefault", desc); } @@ -288,7 +379,7 @@ namespace { } }; // namespace -GsfElectronProducer::GsfElectronProducer(const edm::ParameterSet& cfg, const GsfElectronAlgo::HeavyObjectCache*) +GsfElectronProducer::GsfElectronProducer(const edm::ParameterSet& cfg, const GsfElectronAlgo::HeavyObjectCache* gcache) : cutsCfg_{makeCutsConfiguration(cfg.getParameter("preselection"))}, ecalSeedingParametersChecked_(false), electronPutToken_(produces()), @@ -311,6 +402,19 @@ GsfElectronProducer::GsfElectronProducer(const edm::ParameterSet& cfg, const Gsf if (cfg.getParameter("fillConvVtxFitProb")) inputCfg_.conversions = consumes(cfg.getParameter("conversionsTag")); + // inputs for PFCluster isolation + const edm::ParameterSet& pfECALClusIsolCfg = cfg.getParameter("pfECALClusIsolCfg"); + const edm::ParameterSet& pfHCALClusIsolCfg = cfg.getParameter("pfHCALClusIsolCfg"); + inputCfg_.pfClusterProducer = + consumes(pfECALClusIsolCfg.getParameter("pfClusterProducer")); + inputCfg_.pfClusterProducerHCAL = consumes(pfHCALClusIsolCfg.getParameter("pfClusterProducerHCAL")); + inputCfg_.pfClusterProducerHFEM = consumes(pfHCALClusIsolCfg.getParameter("pfClusterProducerHFEM")); + inputCfg_.pfClusterProducerHFHAD = consumes(pfHCALClusIsolCfg.getParameter("pfClusterProducerHFHAD")); + + // Config for PFID dnn + const auto& pset_dnn = cfg.getParameter("EleDNNPFid"); + dnnPFidEnabled_ = pset_dnn.getParameter("enabled"); + strategyCfg_.useDefaultEnergyCorrection = cfg.getParameter("useDefaultEnergyCorrection"); strategyCfg_.applyPreselection = cfg.getParameter("applyPreselection"); @@ -330,6 +434,7 @@ GsfElectronProducer::GsfElectronProducer(const edm::ParameterSet& cfg, const Gsf strategyCfg_.useEcalRegression = cfg.getParameter("useEcalRegression"); strategyCfg_.useCombinationRegression = cfg.getParameter("useCombinationRegression"); strategyCfg_.fillConvVtxFitProb = cfg.getParameter("fillConvVtxFitProb"); + strategyCfg_.computePfClusterIso = dnnPFidEnabled_; // hcal helpers auto const& psetPreselection = cfg.getParameter("preselection"); @@ -388,6 +493,25 @@ GsfElectronProducer::GsfElectronProducer(const edm::ParameterSet& cfg, const Gsf .vetoClustered = cfg.getParameter("vetoClustered"), .useNumCrystals = cfg.getParameter("useNumCrystals")}; + // isolation + const GsfElectronAlgo::PFClusterIsolationConfiguration pfisoCfg{ + .ecaldrMax = pfECALClusIsolCfg.getParameter("drMax"), + .ecaldrVetoBarrel = pfECALClusIsolCfg.getParameter("drVetoBarrel"), + .ecaldrVetoEndcap = pfECALClusIsolCfg.getParameter("drVetoEndcap"), + .ecaletaStripBarrel = pfECALClusIsolCfg.getParameter("etaStripBarrel"), + .ecaletaStripEndcap = pfECALClusIsolCfg.getParameter("etaStripEndcap"), + .ecalenergyBarrel = pfECALClusIsolCfg.getParameter("energyBarrel"), + .ecalenergyEndcap = pfECALClusIsolCfg.getParameter("energyEndcap"), + .useHF = pfHCALClusIsolCfg.getParameter("useHF"), + .hcaldrMax = pfHCALClusIsolCfg.getParameter("drMax"), + .hcaldrVetoBarrel = pfHCALClusIsolCfg.getParameter("drVetoBarrel"), + .hcaldrVetoEndcap = pfHCALClusIsolCfg.getParameter("drVetoEndcap"), + .hcaletaStripBarrel = pfHCALClusIsolCfg.getParameter("etaStripBarrel"), + .hcaletaStripEndcap = pfHCALClusIsolCfg.getParameter("etaStripEndcap"), + .hcalenergyBarrel = pfHCALClusIsolCfg.getParameter("energyBarrel"), + .hcalenergyEndcap = pfHCALClusIsolCfg.getParameter("energyEndcap"), + .hcaluseEt = pfHCALClusIsolCfg.getParameter("useEt")}; + const RegressionHelper::Configuration regressionCfg{ .ecalRegressionWeightLabels = cfg.getParameter>("ecalRefinedRegressionWeightLabels"), .ecalWeightsFromDB = cfg.getParameter("ecalWeightsFromDB"), @@ -406,6 +530,7 @@ GsfElectronProducer::GsfElectronProducer(const edm::ParameterSet& cfg, const Gsf hcalCfg_, hcalCfgBc_, isoCfg, + pfisoCfg, recHitsCfg, EcalClusterFunctionFactory::get()->create( cfg.getParameter("crackCorrectionFunction"), cfg, consumesCollector()), @@ -415,6 +540,16 @@ GsfElectronProducer::GsfElectronProducer(const edm::ParameterSet& cfg, const Gsf cfg.getParameter("trkIsolHEEP03Cfg"), cfg.getParameter("trkIsolHEEP04Cfg"), consumesCollector()); + + if (dnnPFidEnabled_) { + tfSessions_ = gcache->iElectronDNNEstimator->getSessions(); + } +} + +void GsfElectronProducer::endStream() { + for (auto session : tfSessions_) { + tensorflow::closeSession(session); + } } void GsfElectronProducer::checkEcalSeedingParameters(edm::ParameterSet const& pset) { @@ -586,9 +721,10 @@ void GsfElectronProducer::produce(edm::Event& event, const edm::EventSetup& setu auto electrons = algo_->completeElectrons(event, setup, globalCache()); if (resetMvaValuesUsingPFCandidates_) { const auto gsfMVAInputMap = matchWithPFCandidates(event.get(egmPFCandidateCollection_)); - setMVAOutputs(electrons, globalCache(), event.get(inputCfg_.vtxCollectionTag)); - for (auto& el : electrons) - el.setMvaInput(gsfMVAInputMap.find(el.gsfTrack())->second); // set MVA inputs + for (auto& el : electrons) { + el.setMvaInput(gsfMVAInputMap.find(el.gsfTrack())->second); // set Run2 MVA inputs + } + setMVAOutputs(electrons, globalCache(), event.get(inputCfg_.vtxCollectionTag), dnnPFidEnabled_, tfSessions_); } // all electrons diff --git a/RecoEgamma/EgammaElectronProducers/python/gedGsfElectrons_cfi.py b/RecoEgamma/EgammaElectronProducers/python/gedGsfElectrons_cfi.py index 78fb28793ba9d..9769315525ee6 100644 --- a/RecoEgamma/EgammaElectronProducers/python/gedGsfElectrons_cfi.py +++ b/RecoEgamma/EgammaElectronProducers/python/gedGsfElectrons_cfi.py @@ -19,9 +19,24 @@ "gedelectron_EBUncertainty_offline_v1", "gedelectron_EEUncertainty_offline_v1"], combinationRegressionWeightLabels = ["gedelectron_p4combination_offline"], + + #Activate the evaluation of Egamma PFID DNN + EleDNNPFid= dict( + modelsFiles = [ + "RecoEgamma/ElectronIdentification/data/Ele_PFID_dnn/lowpT/lowpT_modelDNN.pb", + "RecoEgamma/ElectronIdentification/data/Ele_PFID_dnn/highpTEB/highpTEB_modelDNN.pb", + "RecoEgamma/ElectronIdentification/data/Ele_PFID_dnn/highpTEE/highpTEE_modelDNN.pb" + ], + scalersFiles = [ + "RecoEgamma/ElectronIdentification/data/Ele_PFID_dnn/lowpT/lowpT_scaler.txt", + "RecoEgamma/ElectronIdentification/data/Ele_PFID_dnn/highpTEB/highpTEB_scaler.txt", + "RecoEgamma/ElectronIdentification/data/Ele_PFID_dnn/highpTEE/highpTEE_scaler.txt" + ] + ) ) + from Configuration.ProcessModifiers.pp_on_AA_cff import pp_on_AA pp_on_AA.toModify(gedGsfElectronsTmp.preselection, minSCEtBarrel = 15.0) pp_on_AA.toModify(gedGsfElectronsTmp.preselection, minSCEtEndcaps = 15.0) @@ -31,3 +46,10 @@ minSCEtBarrel = 1.0, minSCEtEndcaps = 1.0) egamma_lowPt_exclusive.toModify(gedGsfElectronsTmp, applyPreselection = False) + + +# Activate the Egamma PFID dnn only for Run3 +from Configuration.Eras.Modifier_run3_common_cff import run3_common +run3_common.toModify(gedGsfElectronsTmp.EleDNNPFid, + enabled = True +) \ No newline at end of file diff --git a/RecoEgamma/EgammaElectronProducers/python/gsfElectronProducer_cfi.py b/RecoEgamma/EgammaElectronProducers/python/gsfElectronProducer_cfi.py index 7aa0ae0646090..5a70acfebe75c 100644 --- a/RecoEgamma/EgammaElectronProducers/python/gsfElectronProducer_cfi.py +++ b/RecoEgamma/EgammaElectronProducers/python/gsfElectronProducer_cfi.py @@ -2,10 +2,42 @@ from RecoEgamma.EgammaIsolationAlgos.egammaHBHERecHitThreshold_cff import egammaHBHERecHit +from RecoEgamma.EgammaIsolationAlgos.egammaEcalPFClusterIsolationProducerRecoGsfElectron_cfi import egammaEcalPFClusterIsolationProducerRecoGsfElectron +from RecoEgamma.EgammaIsolationAlgos.egammaHcalPFClusterIsolationProducerRecoGsfElectron_cfi import egammaHcalPFClusterIsolationProducerRecoGsfElectron + import RecoEgamma.EgammaElectronProducers.gsfElectronProducerDefault_cfi as _gsfProd gsfElectronProducer = _gsfProd.gsfElectronProducerDefault.clone( hbheRecHits = egammaHBHERecHit.hbheRecHits, recHitEThresholdHB = egammaHBHERecHit.recHitEThresholdHB, recHitEThresholdHE = egammaHBHERecHit.recHitEThresholdHE, - maxHcalRecHitSeverity = egammaHBHERecHit.maxHcalRecHitSeverity + maxHcalRecHitSeverity = egammaHBHERecHit.maxHcalRecHitSeverity, + + pfECALClusIsolCfg = cms.PSet( + pfClusterProducer = egammaEcalPFClusterIsolationProducerRecoGsfElectron.pfClusterProducer, + drMax = egammaEcalPFClusterIsolationProducerRecoGsfElectron.drMax, + drVetoBarrel = egammaEcalPFClusterIsolationProducerRecoGsfElectron.drVetoBarrel, + drVetoEndcap = egammaEcalPFClusterIsolationProducerRecoGsfElectron.drVetoEndcap, + etaStripBarrel = egammaEcalPFClusterIsolationProducerRecoGsfElectron.etaStripBarrel, + etaStripEndcap = egammaEcalPFClusterIsolationProducerRecoGsfElectron.etaStripEndcap, + energyBarrel = egammaEcalPFClusterIsolationProducerRecoGsfElectron.energyBarrel, + energyEndcap = egammaEcalPFClusterIsolationProducerRecoGsfElectron.energyEndcap + ), + + pfHCALClusIsolCfg = cms.PSet( + + pfClusterProducerHCAL = egammaHcalPFClusterIsolationProducerRecoGsfElectron.pfClusterProducerHCAL, + useHF = egammaHcalPFClusterIsolationProducerRecoGsfElectron.useHF, + pfClusterProducerHFEM = egammaHcalPFClusterIsolationProducerRecoGsfElectron.pfClusterProducerHFEM, + pfClusterProducerHFHAD = egammaHcalPFClusterIsolationProducerRecoGsfElectron.pfClusterProducerHFHAD, + drMax = egammaHcalPFClusterIsolationProducerRecoGsfElectron.drMax, + drVetoBarrel = egammaHcalPFClusterIsolationProducerRecoGsfElectron.drVetoBarrel, + drVetoEndcap = egammaHcalPFClusterIsolationProducerRecoGsfElectron.drVetoEndcap, + etaStripBarrel = egammaHcalPFClusterIsolationProducerRecoGsfElectron.etaStripBarrel, + etaStripEndcap = egammaHcalPFClusterIsolationProducerRecoGsfElectron.etaStripEndcap, + energyBarrel = egammaHcalPFClusterIsolationProducerRecoGsfElectron.energyBarrel, + energyEndcap = egammaHcalPFClusterIsolationProducerRecoGsfElectron.energyEndcap, + useEt = egammaHcalPFClusterIsolationProducerRecoGsfElectron.useEt, + + ) + ) diff --git a/RecoEgamma/EgammaIsolationAlgos/interface/EcalPFClusterIsolation.h b/RecoEgamma/EgammaIsolationAlgos/interface/EcalPFClusterIsolation.h index 8c4d74ebd397f..60ffeff67a0e0 100644 --- a/RecoEgamma/EgammaIsolationAlgos/interface/EcalPFClusterIsolation.h +++ b/RecoEgamma/EgammaIsolationAlgos/interface/EcalPFClusterIsolation.h @@ -28,10 +28,12 @@ class EcalPFClusterIsolation { double energyEndcap); ~EcalPFClusterIsolation(); + double getSum(T1, edm::Handle >); double getSum(T1Ref, edm::Handle >); private: bool computedRVeto(T1Ref candRef, reco::PFClusterRef pfclu); + bool computedRVeto(T1 cand, reco::PFClusterRef pfclu); double drVeto2_; const double drMax_; diff --git a/RecoEgamma/EgammaIsolationAlgos/interface/HcalPFClusterIsolation.h b/RecoEgamma/EgammaIsolationAlgos/interface/HcalPFClusterIsolation.h index e6074880cf171..d41e4b7b49ec5 100644 --- a/RecoEgamma/EgammaIsolationAlgos/interface/HcalPFClusterIsolation.h +++ b/RecoEgamma/EgammaIsolationAlgos/interface/HcalPFClusterIsolation.h @@ -30,6 +30,7 @@ class HcalPFClusterIsolation { ~HcalPFClusterIsolation(); double getSum(const T1Ref candRef, const std::vector>& clusterHandles); + double getSum(const T1 cand, const std::vector>& clusterHandles); private: const double drMax_; diff --git a/RecoEgamma/EgammaIsolationAlgos/src/EcalPFClusterIsolation.cc b/RecoEgamma/EgammaIsolationAlgos/src/EcalPFClusterIsolation.cc index f3115977a8961..2bcae8b02cd7e 100644 --- a/RecoEgamma/EgammaIsolationAlgos/src/EcalPFClusterIsolation.cc +++ b/RecoEgamma/EgammaIsolationAlgos/src/EcalPFClusterIsolation.cc @@ -35,11 +35,40 @@ template EcalPFClusterIsolation::~EcalPFClusterIsolation() {} template -double EcalPFClusterIsolation::getSum(const T1Ref candRef, edm::Handle clusterHandle) { +bool EcalPFClusterIsolation::computedRVeto(T1Ref candRef, reco::PFClusterRef pfclu) { + float dR2 = deltaR2(candRef->eta(), candRef->phi(), pfclu->eta(), pfclu->phi()); + if (dR2 > (drMax_ * drMax_)) + return false; + + if (candRef->superCluster().isNonnull()) { + // Exclude clusters that are part of the candidate + for (reco::CaloCluster_iterator it = candRef->superCluster()->clustersBegin(); + it != candRef->superCluster()->clustersEnd(); + ++it) { + if ((*it)->seed() == pfclu->seed()) { + return false; + } + } + } + + return true; +} + +template <> +bool EcalPFClusterIsolation::computedRVeto(T1Ref candRef, reco::PFClusterRef pfclu) { + float dR2 = deltaR2(candRef->eta(), candRef->phi(), pfclu->eta(), pfclu->phi()); + if (dR2 > (drMax_ * drMax_) || dR2 < drVeto2_) + return false; + else + return true; +} + +template +double EcalPFClusterIsolation::getSum(const T1 cand, edm::Handle clusterHandle) { drVeto2_ = -1.; float etaStrip = -1; - if (fabs(candRef->eta()) < 1.479) { + if (std::abs(cand.eta()) < 1.479) { drVeto2_ = drVetoBarrel_ * drVetoBarrel_; etaStrip = etaStripBarrel_; } else { @@ -51,18 +80,18 @@ double EcalPFClusterIsolation::getSum(const T1Ref candRef, edm::Handlesize(); i++) { reco::PFClusterRef pfclu(clusterHandle, i); - if (fabs(candRef->eta()) < 1.479) { - if (fabs(pfclu->pt()) < energyBarrel_) + if (std::abs(cand.eta()) < 1.479) { + if (std::abs(pfclu->pt()) < energyBarrel_) continue; } else { - if (fabs(pfclu->energy()) < energyEndcap_) + if (std::abs(pfclu->energy()) < energyEndcap_) continue; } - float dEta = fabs(candRef->eta() - pfclu->eta()); + float dEta = std::abs(cand.eta() - pfclu->eta()); if (dEta < etaStrip) continue; - if (not computedRVeto(candRef, pfclu)) + if (not computedRVeto(cand, pfclu)) continue; etSum += pfclu->pt(); @@ -72,15 +101,19 @@ double EcalPFClusterIsolation::getSum(const T1Ref candRef, edm::Handle -bool EcalPFClusterIsolation::computedRVeto(T1Ref candRef, reco::PFClusterRef pfclu) { - float dR2 = deltaR2(candRef->eta(), candRef->phi(), pfclu->eta(), pfclu->phi()); +double EcalPFClusterIsolation::getSum(T1Ref ref, edm::Handle > clusts) { + return getSum(*ref, clusts); +} + +template +bool EcalPFClusterIsolation::computedRVeto(T1 cand, reco::PFClusterRef pfclu) { + float dR2 = deltaR2(cand.eta(), cand.phi(), pfclu->eta(), pfclu->phi()); if (dR2 > (drMax_ * drMax_)) return false; - if (candRef->superCluster().isNonnull()) { + if (cand.superCluster().isNonnull()) { // Exclude clusters that are part of the candidate - for (reco::CaloCluster_iterator it = candRef->superCluster()->clustersBegin(); - it != candRef->superCluster()->clustersEnd(); + for (reco::CaloCluster_iterator it = cand.superCluster()->clustersBegin(); it != cand.superCluster()->clustersEnd(); ++it) { if ((*it)->seed() == pfclu->seed()) { return false; @@ -91,15 +124,6 @@ bool EcalPFClusterIsolation::computedRVeto(T1Ref candRef, reco::PFClusterRef return true; } -template <> -bool EcalPFClusterIsolation::computedRVeto(T1Ref candRef, reco::PFClusterRef pfclu) { - float dR2 = deltaR2(candRef->eta(), candRef->phi(), pfclu->eta(), pfclu->phi()); - if (dR2 > (drMax_ * drMax_) || dR2 < drVeto2_) - return false; - else - return true; -} - template class EcalPFClusterIsolation; template class EcalPFClusterIsolation; template class EcalPFClusterIsolation; diff --git a/RecoEgamma/EgammaIsolationAlgos/src/HcalPFClusterIsolation.cc b/RecoEgamma/EgammaIsolationAlgos/src/HcalPFClusterIsolation.cc index dacf23a75dd38..b4d03e5a0a366 100644 --- a/RecoEgamma/EgammaIsolationAlgos/src/HcalPFClusterIsolation.cc +++ b/RecoEgamma/EgammaIsolationAlgos/src/HcalPFClusterIsolation.cc @@ -29,13 +29,14 @@ template HcalPFClusterIsolation::~HcalPFClusterIsolation() {} template -double HcalPFClusterIsolation::getSum(const T1Ref candRef, +double HcalPFClusterIsolation::getSum(const T1 cand, const std::vector>& clusterHandles) { double etSum = 0.; + double candAbsEta = std::abs(cand.eta()); float etaStrip = 0; float dRVeto = 0; - if (fabs(candRef->eta()) < 1.479) { + if (candAbsEta < 1.479) { dRVeto = drVetoBarrel_; etaStrip = etaStripBarrel_; } else { @@ -47,20 +48,20 @@ double HcalPFClusterIsolation::getSum(const T1Ref candRef, for (unsigned i = 0; i < clusterHandles[nHandle]->size(); i++) { const reco::PFClusterRef pfclu(clusterHandles[nHandle], i); - if (fabs(candRef->eta()) < 1.479) { - if (fabs(pfclu->pt()) < energyBarrel_) + if (candAbsEta < 1.479) { + if (std::abs(pfclu->pt()) < energyBarrel_) continue; } else { - if (fabs(pfclu->energy()) < energyEndcap_) + if (std::abs(pfclu->energy()) < energyEndcap_) continue; } - float dEta = fabs(candRef->eta() - pfclu->eta()); + float dEta = std::abs(cand.eta() - pfclu->eta()); if (dEta < etaStrip) continue; - float dR = deltaR(candRef->eta(), candRef->phi(), pfclu->eta(), pfclu->phi()); - if (dR > drMax_ || dR < dRVeto) + float dR2 = deltaR2(cand.eta(), cand.phi(), pfclu->eta(), pfclu->phi()); + if (dR2 > (drMax_ * drMax_) || dR2 < (dRVeto * dRVeto)) continue; if (useEt_) @@ -73,6 +74,12 @@ double HcalPFClusterIsolation::getSum(const T1Ref candRef, return etSum; } +template +double HcalPFClusterIsolation::getSum(T1Ref ref, + const std::vector>& clusterHandles) { + return getSum(*ref, clusterHandles); +} + template class HcalPFClusterIsolation; template class HcalPFClusterIsolation; template class HcalPFClusterIsolation; diff --git a/RecoEgamma/EgammaPhotonProducers/BuildFile.xml b/RecoEgamma/EgammaPhotonProducers/BuildFile.xml index d63630dd1b231..ef90a45005c96 100644 --- a/RecoEgamma/EgammaPhotonProducers/BuildFile.xml +++ b/RecoEgamma/EgammaPhotonProducers/BuildFile.xml @@ -22,4 +22,5 @@ + diff --git a/RecoEgamma/EgammaPhotonProducers/python/gedPhotonSequence_cff.py b/RecoEgamma/EgammaPhotonProducers/python/gedPhotonSequence_cff.py index eee0196b55b16..0718d38d2b983 100644 --- a/RecoEgamma/EgammaPhotonProducers/python/gedPhotonSequence_cff.py +++ b/RecoEgamma/EgammaPhotonProducers/python/gedPhotonSequence_cff.py @@ -13,7 +13,14 @@ photonProducer = "gedPhotonCore", candidateP4type = "fromEcalEnergy", outputPhotonCollection = "", - reconstructionStep = "tmp" + reconstructionStep = "tmp", + PhotonDNNPFid = dict( + modelsFiles = [ "RecoEgamma/PhotonIdentification/data/Photon_PFID_dnn/EB/EB_modelDNN.pb", + "RecoEgamma/PhotonIdentification/data/Photon_PFID_dnn/EE/EE_modelDNN.pb"], + scalersFiles = [ + "RecoEgamma/PhotonIdentification/data/Photon_PFID_dnn/EB/EB_scaler.txt", + "RecoEgamma/PhotonIdentification/data/Photon_PFID_dnn/EE/EE_scaler.txt"] + ) ) del gedPhotonsTmp.regressionConfig @@ -27,13 +34,13 @@ pfECALClusIsolation = cms.InputTag("photonEcalPFClusterIsolationProducer"), pfHCALClusIsolation = cms.InputTag("photonHcalPFClusterIsolationProducer"), pfIsolCfg = cms.PSet( - chargedHadronIso = cms.InputTag("photonIDValueMaps","phoChargedIsolation"), - neutralHadronIso = cms.InputTag("photonIDValueMaps","phoNeutralHadronIsolation"), - photonIso = cms.InputTag("photonIDValueMaps","phoPhotonIsolation"), - chargedHadronWorstVtxIso = cms.InputTag("photonIDValueMaps","phoWorstChargedIsolation"), - chargedHadronWorstVtxGeomVetoIso = cms.InputTag("photonIDValueMaps","phoWorstChargedIsolationConeVeto"), - chargedHadronPFPVIso = cms.InputTag("egmPhotonIsolationCITK:h+-DR030-"), - ) + chargedHadronIso = cms.InputTag("photonIDValueMaps","phoChargedIsolation"), + neutralHadronIso = cms.InputTag("photonIDValueMaps","phoNeutralHadronIsolation"), + photonIso = cms.InputTag("photonIDValueMaps","phoPhotonIsolation"), + chargedHadronWorstVtxIso = cms.InputTag("photonIDValueMaps","phoWorstChargedIsolation"), + chargedHadronWorstVtxGeomVetoIso = cms.InputTag("photonIDValueMaps","phoWorstChargedIsolationConeVeto"), + chargedHadronPFPVIso = cms.InputTag("egmPhotonIsolationCITK:h+-DR030-"), + ), ) gedPhotonTask = cms.Task(gedPhotons) gedPhotonSequence = cms.Sequence(gedPhotonTask) @@ -45,3 +52,10 @@ egamma_lowPt_exclusive.toModify(gedPhotonsTmp, minSCEtBarrel = 1.0, minSCEtEndcap = 1.0) + + +# Activate the Egamma PFID dnn only for Run3 +from Configuration.Eras.Modifier_run3_common_cff import run3_common +run3_common.toModify(gedPhotonsTmp.PhotonDNNPFid, + enabled = True +) \ No newline at end of file diff --git a/RecoEgamma/EgammaPhotonProducers/python/gedPhotons_cfi.py b/RecoEgamma/EgammaPhotonProducers/python/gedPhotons_cfi.py index 65685cd1a3473..4b89499cb2aec 100644 --- a/RecoEgamma/EgammaPhotonProducers/python/gedPhotons_cfi.py +++ b/RecoEgamma/EgammaPhotonProducers/python/gedPhotons_cfi.py @@ -9,6 +9,10 @@ from RecoEgamma.EgammaTools.regressionModifier_cfi import * from RecoEgamma.EgammaIsolationAlgos.egammaHBHERecHitThreshold_cff import egammaHBHERecHit + +from RecoEgamma.EgammaIsolationAlgos.egammaEcalPFClusterIsolationProducerRecoPhoton_cfi import egammaEcalPFClusterIsolationProducerRecoPhoton +from RecoEgamma.EgammaIsolationAlgos.egammaHcalPFClusterIsolationProducerRecoPhoton_cfi import egammaHcalPFClusterIsolationProducerRecoPhoton + # # producer for photons # @@ -94,7 +98,50 @@ RecHitSeverityToBeExcludedEB = cleanedHybridSuperClusters.RecHitSeverityToBeExcluded, RecHitFlagToBeExcludedEE = multi5x5BasicClustersCleaned.RecHitFlagToBeExcluded, RecHitSeverityToBeExcludedEE = cleanedHybridSuperClusters.RecHitSeverityToBeExcluded, - checkHcalStatus = cms.bool(True) + checkHcalStatus = cms.bool(True), + PhotonDNNPFid = cms.PSet( + enabled = cms.bool(False), + inputTensorName = cms.string("FirstLayer_input"), + outputTensorName = cms.string("sequential/FinalLayer/Sigmoid"), + modelsFiles = cms.vstring( + 'RecoEgamma/PhotonIdentification/data/Photon_PFID_dnn/EB/EB_modelDNN.pb', + 'RecoEgamma/PhotonIdentification/data/Photon_PFID_dnn/EE/EE_modelDNN.pb'), + scalersFiles = cms.vstring( + 'RecoEgamma/PhotonIdentification/data/Photon_PFID_dnn/EB/EB_scaler.txt', + 'RecoEgamma/PhotonIdentification/data/Photon_PFID_dnn/EE/EE_scaler.txt' + ), + outputDim = cms.uint32(1), + useEBModelInGap = cms.bool(True) + ), + + pfECALClusIsolCfg = cms.PSet( + pfClusterProducer = egammaEcalPFClusterIsolationProducerRecoPhoton.pfClusterProducer, + drMax = egammaEcalPFClusterIsolationProducerRecoPhoton.drMax, + drVetoBarrel = egammaEcalPFClusterIsolationProducerRecoPhoton.drVetoBarrel, + drVetoEndcap = egammaEcalPFClusterIsolationProducerRecoPhoton.drVetoEndcap, + etaStripBarrel = egammaEcalPFClusterIsolationProducerRecoPhoton.etaStripBarrel, + etaStripEndcap = egammaEcalPFClusterIsolationProducerRecoPhoton.etaStripEndcap, + energyBarrel = egammaEcalPFClusterIsolationProducerRecoPhoton.energyBarrel, + energyEndcap = egammaEcalPFClusterIsolationProducerRecoPhoton.energyEndcap + ), + + pfHCALClusIsolCfg = cms.PSet( + + pfClusterProducerHCAL = egammaHcalPFClusterIsolationProducerRecoPhoton.pfClusterProducerHCAL, + useHF = egammaHcalPFClusterIsolationProducerRecoPhoton.useHF, + pfClusterProducerHFEM = egammaHcalPFClusterIsolationProducerRecoPhoton.pfClusterProducerHFEM, + pfClusterProducerHFHAD = egammaHcalPFClusterIsolationProducerRecoPhoton.pfClusterProducerHFHAD, + drMax = egammaHcalPFClusterIsolationProducerRecoPhoton.drMax, + drVetoBarrel = egammaHcalPFClusterIsolationProducerRecoPhoton.drVetoBarrel, + drVetoEndcap = egammaHcalPFClusterIsolationProducerRecoPhoton.drVetoEndcap, + etaStripBarrel = egammaHcalPFClusterIsolationProducerRecoPhoton.etaStripBarrel, + etaStripEndcap = egammaHcalPFClusterIsolationProducerRecoPhoton.etaStripEndcap, + energyBarrel = egammaHcalPFClusterIsolationProducerRecoPhoton.energyBarrel, + energyEndcap = egammaHcalPFClusterIsolationProducerRecoPhoton.energyEndcap, + useEt = egammaHcalPFClusterIsolationProducerRecoPhoton.useEt, + + ) + ) diff --git a/RecoEgamma/EgammaPhotonProducers/python/photons_cfi.py b/RecoEgamma/EgammaPhotonProducers/python/photons_cfi.py index 85b4699c819c9..7b8e9c69bbc35 100644 --- a/RecoEgamma/EgammaPhotonProducers/python/photons_cfi.py +++ b/RecoEgamma/EgammaPhotonProducers/python/photons_cfi.py @@ -5,6 +5,10 @@ from RecoEcal.EgammaClusterProducers.hybridSuperClusters_cfi import * from RecoEcal.EgammaClusterProducers.multi5x5BasicClusters_cfi import * from RecoEgamma.EgammaIsolationAlgos.egammaHBHERecHitThreshold_cff import egammaHBHERecHit + +from RecoEgamma.EgammaIsolationAlgos.egammaEcalPFClusterIsolationProducerRecoPhoton_cfi import egammaEcalPFClusterIsolationProducerRecoPhoton +from RecoEgamma.EgammaIsolationAlgos.egammaHcalPFClusterIsolationProducerRecoPhoton_cfi import egammaHcalPFClusterIsolationProducerRecoPhoton + # # producer for photons # @@ -85,7 +89,49 @@ RecHitSeverityToBeExcludedEB = cleanedHybridSuperClusters.RecHitSeverityToBeExcluded, RecHitFlagToBeExcludedEE = multi5x5BasicClustersCleaned.RecHitFlagToBeExcluded, RecHitSeverityToBeExcludedEE = cleanedHybridSuperClusters.RecHitSeverityToBeExcluded, - checkHcalStatus = cms.bool(True) + checkHcalStatus = cms.bool(True), + PhotonDNNPFid = cms.PSet( + enabled = cms.bool(False), + inputTensorName = cms.string("FirstLayer_input"), + outputTensorName = cms.string("sequential/LastLayer/Softmax"), + modelsFiles = cms.vstring( + 'RecoEgamma/PhotonIdentification/data/Photon_PFID_dnn/model_barrel.pb', + 'RecoEgamma/PhotonIdentification/data/Photon_PFID_dnn/model_endcap.pb'), + scalersFiles = cms.vstring( + 'RecoEgamma/PhotonIdentification/data/Photon_PFID_dnn/PhotonDNNScaler.txt', + 'RecoEgamma/PhotonIdentification/data/Photon_PFID_dnn/PhotonDNNScaler.txt' + ), + outputDim = cms.uint32(1), + useEBModelInGap = cms.bool(True) + ), + pfECALClusIsolCfg = cms.PSet( + pfClusterProducer = egammaEcalPFClusterIsolationProducerRecoPhoton.pfClusterProducer, + drMax = egammaEcalPFClusterIsolationProducerRecoPhoton.drMax, + drVetoBarrel = egammaEcalPFClusterIsolationProducerRecoPhoton.drVetoBarrel, + drVetoEndcap = egammaEcalPFClusterIsolationProducerRecoPhoton.drVetoEndcap, + etaStripBarrel = egammaEcalPFClusterIsolationProducerRecoPhoton.etaStripBarrel, + etaStripEndcap = egammaEcalPFClusterIsolationProducerRecoPhoton.etaStripEndcap, + energyBarrel = egammaEcalPFClusterIsolationProducerRecoPhoton.energyBarrel, + energyEndcap = egammaEcalPFClusterIsolationProducerRecoPhoton.energyEndcap + ), + + pfHCALClusIsolCfg = cms.PSet( + + pfClusterProducerHCAL = egammaHcalPFClusterIsolationProducerRecoPhoton.pfClusterProducerHCAL, + useHF = egammaHcalPFClusterIsolationProducerRecoPhoton.useHF, + pfClusterProducerHFEM = egammaHcalPFClusterIsolationProducerRecoPhoton.pfClusterProducerHFEM, + pfClusterProducerHFHAD = egammaHcalPFClusterIsolationProducerRecoPhoton.pfClusterProducerHFHAD, + drMax = egammaHcalPFClusterIsolationProducerRecoPhoton.drMax, + drVetoBarrel = egammaHcalPFClusterIsolationProducerRecoPhoton.drVetoBarrel, + drVetoEndcap = egammaHcalPFClusterIsolationProducerRecoPhoton.drVetoEndcap, + etaStripBarrel = egammaHcalPFClusterIsolationProducerRecoPhoton.etaStripBarrel, + etaStripEndcap = egammaHcalPFClusterIsolationProducerRecoPhoton.etaStripEndcap, + energyBarrel = egammaHcalPFClusterIsolationProducerRecoPhoton.energyBarrel, + energyEndcap = egammaHcalPFClusterIsolationProducerRecoPhoton.energyEndcap, + useEt = egammaHcalPFClusterIsolationProducerRecoPhoton.useEt, + + ) + ) photonsHGC = photons.clone( diff --git a/RecoEgamma/EgammaPhotonProducers/src/GEDPhotonProducer.cc b/RecoEgamma/EgammaPhotonProducers/src/GEDPhotonProducer.cc index 6f953ca768599..337cea3032416 100644 --- a/RecoEgamma/EgammaPhotonProducers/src/GEDPhotonProducer.cc +++ b/RecoEgamma/EgammaPhotonProducers/src/GEDPhotonProducer.cc @@ -46,13 +46,42 @@ #include "CondFormats/DataRecord/interface/EcalPFRecHitThresholdsRcd.h" #include "RecoEcal/EgammaCoreTools/interface/EgammaLocalCovParamDefaults.h" #include "RecoEgamma/EgammaElectronAlgos/interface/ElectronHcalHelper.h" +#include "RecoEgamma/PhotonIdentification/interface/PhotonDNNEstimator.h" +#include "PhysicsTools/TensorFlow/interface/TensorFlow.h" +#include "RecoEgamma/EgammaIsolationAlgos/interface/EcalPFClusterIsolation.h" +#include "RecoEgamma/EgammaIsolationAlgos/interface/HcalPFClusterIsolation.h" -class GEDPhotonProducer : public edm::stream::EDProducer<> { +class CacheData { public: - GEDPhotonProducer(const edm::ParameterSet& ps); + CacheData(const edm::ParameterSet& conf) { + // Here we will have to load the DNN PFID if present in the config + egammaTools::DNNConfiguration config; + const auto& pset_dnn = conf.getParameter("PhotonDNNPFid"); + const auto dnnEnabled = pset_dnn.getParameter("enabled"); + if (dnnEnabled) { + config.inputTensorName = pset_dnn.getParameter("inputTensorName"); + config.outputTensorName = pset_dnn.getParameter("outputTensorName"); + config.modelsFiles = pset_dnn.getParameter>("modelsFiles"); + config.scalersFiles = pset_dnn.getParameter>("scalersFiles"); + config.outputDim = pset_dnn.getParameter("outputDim"); + const auto useEBModelInGap = pset_dnn.getParameter("useEBModelInGap"); + photonDNNEstimator = std::make_unique(config, useEBModelInGap); + } + } + std::unique_ptr photonDNNEstimator; +}; + +class GEDPhotonProducer : public edm::stream::EDProducer> { +public: + GEDPhotonProducer(const edm::ParameterSet& ps, const CacheData* gcache); void produce(edm::Event& evt, const edm::EventSetup& es) override; + static std::unique_ptr initializeGlobalCache(const edm::ParameterSet&); + static void globalEndJob(const CacheData*){}; + + void endStream() override; + private: class RecoStepInfo { public: @@ -169,6 +198,38 @@ class GEDPhotonProducer : public edm::stream::EDProducer<> { std::unique_ptr hcalHelperCone_; std::unique_ptr hcalHelperBc_; bool hcalRun2EffDepth_; + + // DNN for PFID photon enabled + bool dnnPFidEnabled_; + std::vector tfSessions_; + + double ecaldrMax_; + double ecaldrVetoBarrel_; + double ecaldrVetoEndcap_; + double ecaletaStripBarrel_; + double ecaletaStripEndcap_; + double ecalenergyBarrel_; + double ecalenergyEndcap_; + typedef EcalPFClusterIsolation PhotonEcalPFClusterIsolation; + std::unique_ptr ecalisoAlgo = nullptr; + edm::EDGetTokenT pfClusterProducer_; + + bool useHF_; + double hcaldrMax_; + double hcaldrVetoBarrel_; + double hcaldrVetoEndcap_; + double hcaletaStripBarrel_; + double hcaletaStripEndcap_; + double hcalenergyBarrel_; + double hcalenergyEndcap_; + double hcaluseEt_; + + edm::EDGetTokenT pfClusterProducerHCAL_; + edm::EDGetTokenT pfClusterProducerHFEM_; + edm::EDGetTokenT pfClusterProducerHFHAD_; + + typedef HcalPFClusterIsolation PhotonHcalPFClusterIsolation; + std::unique_ptr hcalisoAlgo = nullptr; }; #include "FWCore/Framework/interface/MakerMacros.h" @@ -196,7 +257,7 @@ GEDPhotonProducer::RecoStepInfo::RecoStepInfo(const std::string& step) : flags_( } } -GEDPhotonProducer::GEDPhotonProducer(const edm::ParameterSet& config) +GEDPhotonProducer::GEDPhotonProducer(const edm::ParameterSet& config, const CacheData* gcache) : photonProducer_{config.getParameter("photonProducer")}, ecalClusterESGetTokens_{consumesCollector()}, recoStep_(config.getParameter("reconstructionStep")), @@ -228,6 +289,7 @@ GEDPhotonProducer::GEDPhotonProducer(const edm::ParameterSet& config) if (config.exists("pfHCALClusIsolation")) { phoPFHCALClusIsolationToken_ = consumes(config.getParameter("pfHCALClusIsolation")); } + } else { photonCoreProducerT_ = consumes(photonProducer_); } @@ -308,8 +370,6 @@ GEDPhotonProducer::GEDPhotonProducer(const edm::ParameterSet& config) hcalRun2EffDepth_ = config.getParameter("hcalRun2EffDepth"); - //AA - // cut values for pre-selection preselCutValuesBarrel_ = {config.getParameter("minSCEtBarrel"), config.getParameter("maxHoverEBarrel"), @@ -351,11 +411,54 @@ GEDPhotonProducer::GEDPhotonProducer(const edm::ParameterSet& config) photonMIPHaloTagger_->setup(mipVariableSet, consumesCollector()); } + ///Get the set for PF cluster isolation calculator + const edm::ParameterSet& pfECALClusIsolCfg = config.getParameter("pfECALClusIsolCfg"); + pfClusterProducer_ = + consumes(pfECALClusIsolCfg.getParameter("pfClusterProducer")); + ecaldrMax_ = pfECALClusIsolCfg.getParameter("drMax"); + ecaldrVetoBarrel_ = pfECALClusIsolCfg.getParameter("drVetoBarrel"); + ecaldrVetoEndcap_ = pfECALClusIsolCfg.getParameter("drVetoEndcap"); + ecaletaStripBarrel_ = pfECALClusIsolCfg.getParameter("etaStripBarrel"); + ecaletaStripEndcap_ = pfECALClusIsolCfg.getParameter("etaStripEndcap"); + ecalenergyBarrel_ = pfECALClusIsolCfg.getParameter("energyBarrel"); + ecalenergyEndcap_ = pfECALClusIsolCfg.getParameter("energyEndcap"); + + const edm::ParameterSet& pfHCALClusIsolCfg = config.getParameter("pfHCALClusIsolCfg"); + pfClusterProducerHCAL_ = consumes(pfHCALClusIsolCfg.getParameter("pfClusterProducerHCAL")); + pfClusterProducerHFEM_ = consumes(pfHCALClusIsolCfg.getParameter("pfClusterProducerHFEM")); + pfClusterProducerHFHAD_ = consumes(pfHCALClusIsolCfg.getParameter("pfClusterProducerHFHAD")); + useHF_ = pfHCALClusIsolCfg.getParameter("useHF"); + hcaldrMax_ = pfHCALClusIsolCfg.getParameter("drMax"); + hcaldrVetoBarrel_ = pfHCALClusIsolCfg.getParameter("drVetoBarrel"); + hcaldrVetoEndcap_ = pfHCALClusIsolCfg.getParameter("drVetoEndcap"); + hcaletaStripBarrel_ = pfHCALClusIsolCfg.getParameter("etaStripBarrel"); + hcaletaStripEndcap_ = pfHCALClusIsolCfg.getParameter("etaStripEndcap"); + hcalenergyBarrel_ = pfHCALClusIsolCfg.getParameter("energyBarrel"); + hcalenergyEndcap_ = pfHCALClusIsolCfg.getParameter("energyEndcap"); + hcaluseEt_ = pfHCALClusIsolCfg.getParameter("useEt"); + // Register the product produces(photonCollection_); if (not pfEgammaCandidates_.isUninitialized()) { produces>(valueMapPFCandPhoton_); } + + const auto& pset_dnn = config.getParameter("PhotonDNNPFid"); + dnnPFidEnabled_ = pset_dnn.getParameter("enabled"); + if (dnnPFidEnabled_) { + tfSessions_ = gcache->photonDNNEstimator->getSessions(); + } +} + +std::unique_ptr GEDPhotonProducer::initializeGlobalCache(const edm::ParameterSet& config) { + // this method is supposed to create, initialize and return a CacheData instance + return std::make_unique(config); +} + +void GEDPhotonProducer::endStream() { + for (auto session : tfSessions_) { + tensorflow::closeSession(session); + } } void GEDPhotonProducer::produce(edm::Event& theEvent, const edm::EventSetup& eventSetup) { @@ -464,6 +567,24 @@ void GEDPhotonProducer::produce(edm::Event& theEvent, const edm::EventSetup& eve photonEnergyCorrector_->gedRegression()->setEventContent(eventSetup); } + ///PF ECAL cluster based isolations + ecalisoAlgo = std::make_unique(ecaldrMax_, + ecaldrVetoBarrel_, + ecaldrVetoEndcap_, + ecaletaStripBarrel_, + ecaletaStripEndcap_, + ecalenergyBarrel_, + ecalenergyEndcap_); + + hcalisoAlgo = std::make_unique(hcaldrMax_, + hcaldrVetoBarrel_, + hcaldrVetoEndcap_, + hcaletaStripBarrel_, + hcaletaStripEndcap_, + hcalenergyBarrel_, + hcalenergyEndcap_, + hcaluseEt_); + int iSC = 0; // index in photon collection // Loop over barrel and endcap SC collections and fill the photon collection if (validPhotonCoreHandle) @@ -631,6 +752,7 @@ void GEDPhotonProducer::fillPhotonCollection(edm::Event& evt, float full5x5_sigmaEtaEta = sqrt(full5x5_cov[0]); float full5x5_sigmaIetaIeta = sqrt(full5x5_locCov[0]); + float full5x5_sigmaIetaIphi = full5x5_locCov[1]; // compute position of ECAL shower math::XYZPoint caloPosition = scRef->position(); @@ -673,6 +795,7 @@ void GEDPhotonProducer::fillPhotonCollection(edm::Event& evt, for (uint id = 0; id < showerShape.hcalOverEcal.size(); ++id) { showerShape.hcalOverEcal[id] = (hcalHelperCone != nullptr) ? hcalHelperCone->hcalESum(*scRef, id + 1) / scRef->energy() : 0.f; + showerShape.hcalOverEcalBc[id] = (hcalHelperBc != nullptr) ? hcalHelperBc->hcalESum(*scRef, id + 1) / scRef->energy() : 0.f; } @@ -749,7 +872,7 @@ void GEDPhotonProducer::fillPhotonCollection(edm::Event& evt, full5x5_showerShape.sigmaIetaIeta = full5x5_sigmaIetaIeta; /// fill extra full5x5 shower shapes const float full5x5_spp = (!edm::isFinite(full5x5_locCov[2]) ? 0. : std::sqrt(full5x5_locCov[2])); - const float full5x5_sep = full5x5_locCov[1]; + const float full5x5_sep = full5x5_sigmaIetaIphi; full5x5_showerShape.sigmaIetaIphi = full5x5_sep; full5x5_showerShape.sigmaIphiIphi = full5x5_spp; full5x5_showerShape.e2nd = (hits != nullptr ? noZS::EcalClusterTools::e2nd(*(scRef->seed()), hits) : 0.f); @@ -793,6 +916,27 @@ void GEDPhotonProducer::fillPhotonCollection(edm::Event& evt, full5x5_showerShape.pre7DepthHcal = false; newCandidate.full5x5_setShowerShapeVariables(full5x5_showerShape); + //get the pointer for the photon object + edm::Ptr photonPtr(photonCoreHandle, lSC); + + // New in CMSSW_12_1_0 for PFID with DNNs + // The PFIso values are computed in the first loop on gedPhotonsTmp to make them available as DNN inputs. + // They are computed with the same inputs and algo as the final PFiso variables computed in the second loop after PF. + // Get PFClusters for PFID only if the PFID DNN evaluation is enabled + if (dnnPFidEnabled_) { + auto clusterHandle = evt.getHandle(pfClusterProducer_); + std::vector> clusterHandles{evt.getHandle(pfClusterProducerHCAL_)}; + if (useHF_) { + clusterHandles.push_back(evt.getHandle(pfClusterProducerHFEM_)); + clusterHandles.push_back(evt.getHandle(pfClusterProducerHFHAD_)); + } + reco::Photon::PflowIsolationVariables pfIso; + pfIso.sumEcalClusterEt = ecalisoAlgo->getSum(newCandidate, clusterHandle); + pfIso.sumHcalClusterEt = hcalisoAlgo->getSum(newCandidate, clusterHandles); + + newCandidate.setPflowIsolationVariables(pfIso); + } + /// get ecal photon specific corrected energy /// plus values from regressions and store them in the Photon // Photon candidate takes by default (set in photons_cfi.py) @@ -854,6 +998,20 @@ void GEDPhotonProducer::fillPhotonCollection(edm::Event& evt, if (isLooseEM) outputPhotonCollection.push_back(newCandidate); } + + if (dnnPFidEnabled_) { + // Here send the list of photons to the PhotonDNNEstimator and get back the values for all the photons in one go + LogDebug("GEDPhotonProducer") << "Getting DNN PFId for photons"; + const auto& dnn_photon_pfid = globalCache()->photonDNNEstimator->evaluate(outputPhotonCollection, tfSessions_); + size_t ipho = 0; + for (auto& photon : outputPhotonCollection) { + const auto& values = dnn_photon_pfid[ipho]; + reco::Photon::PflowIDVariables pfID; + pfID.dnn = values[0]; + photon.setPflowIDVariables(pfID); + ipho++; + } + } } void GEDPhotonProducer::fillPhotonCollection(edm::Event& evt, @@ -895,12 +1053,13 @@ void GEDPhotonProducer::fillPhotonCollection(edm::Event& evt, if (parentSCRef.isNonnull() && ptFast(parentSCRef->energy(), parentSCRef->position(), {0, 0, 0}) <= preselCutValues[0]) continue; + reco::Photon newCandidate(*phoRef); iSC++; - // Calculate the PF isolation and ID - for the time being there is no calculation. Only the setting + // Calculate the PF isolation reco::Photon::PflowIsolationVariables pfIso; - reco::Photon::PflowIDVariables pfID; + // The PFID are not recomputed since they have been already computed in the first loop with the DNN //get the pointer for the photon object edm::Ptr photonPtr(photonHandle, lSC); @@ -915,11 +1074,10 @@ void GEDPhotonProducer::fillPhotonCollection(edm::Event& evt, } //OOT photons in legacy 80X reminiAOD workflow dont have pf cluster isolation embeded into them at this stage + // They have been already computed in the first loop on gedPhotonsTmp but better to compute them again here. pfIso.sumEcalClusterEt = !phoPFECALClusIsolationToken_.isUninitialized() ? (*pfEcalClusters)[photonPtr] : 0.; pfIso.sumHcalClusterEt = !phoPFHCALClusIsolationToken_.isUninitialized() ? (*pfHcalClusters)[photonPtr] : 0.; - newCandidate.setPflowIsolationVariables(pfIso); - newCandidate.setPflowIDVariables(pfID); // do the regression photonEnergyCorrector_->calculate(evt, newCandidate, subdet, vertexCollection, es); diff --git a/RecoEgamma/EgammaTools/BuildFile.xml b/RecoEgamma/EgammaTools/BuildFile.xml index 6a5f6b7815913..24be904030635 100644 --- a/RecoEgamma/EgammaTools/BuildFile.xml +++ b/RecoEgamma/EgammaTools/BuildFile.xml @@ -17,6 +17,7 @@ + diff --git a/RecoEgamma/EgammaTools/interface/EgammaDNNHelper.h b/RecoEgamma/EgammaTools/interface/EgammaDNNHelper.h new file mode 100644 index 0000000000000..2752d3475bce3 --- /dev/null +++ b/RecoEgamma/EgammaTools/interface/EgammaDNNHelper.h @@ -0,0 +1,74 @@ +#ifndef RecoEgamma_ElectronTools_EgammaDNNHelper_h +#define RecoEgamma_ElectronTools_EgammaDNNHelper_h + +#include "PhysicsTools/TensorFlow/interface/TensorFlow.h" +#include +#include +#include +#include + +//author: Davide Valsecchi +//description: +// Handles Tensorflow DNN graphs and variables scaler configuration. +// To be used for PFID egamma DNNs + +namespace egammaTools { + + struct DNNConfiguration { + std::string inputTensorName; + std::string outputTensorName; + std::vector modelsFiles; + std::vector scalersFiles; + uint outputDim = 1; + }; + + struct ScalerConfiguration { + /* Each input variable is represented by the tuple + * The standardization_type can be: + * 0 = Do not scale the variable + * 1 = standard norm. par1=mean, par2=std + * 2 = MinMax. par1=min, par2=max */ + std::string varName; + uint type; + float par1; + float par2; + }; + + // Model for function to be used on the specific candidate to get the model + // index to be used for the evaluation. + typedef std::function&)> ModelSelector; + + class EgammaDNNHelper { + public: + EgammaDNNHelper(const DNNConfiguration&, const ModelSelector& sel, const std::vector& availableVars); + + std::vector getSessions() const; + // Function getting the input vector for a specific electron, already scaled + // together with the model index it has to be used. + // The model index is determined by the ModelSelector functor passed in the constructor + // which has access to all the variables. + std::pair> getScaledInputs(const std::map& variables) const; + + std::vector> evaluate(const std::vector>& candidates, + const std::vector& sessions) const; + + private: + void initTensorFlowGraphs(); + void initScalerFiles(const std::vector& availableVars); + + const DNNConfiguration cfg_; + const ModelSelector modelSelector_; + // Number of models handled by the object + uint nModels_; + // Number of inputs for each loaded model + std::vector nInputs_; + + std::vector> graphDefs_; + + // List of input variables for each of the model; + std::vector> featuresMap_; + }; + +}; // namespace egammaTools + +#endif diff --git a/RecoEgamma/EgammaTools/src/EgammaDNNHelper.cc b/RecoEgamma/EgammaTools/src/EgammaDNNHelper.cc new file mode 100644 index 0000000000000..f2cd0a5511c9d --- /dev/null +++ b/RecoEgamma/EgammaTools/src/EgammaDNNHelper.cc @@ -0,0 +1,182 @@ + +#include "RecoEgamma/EgammaTools/interface/EgammaDNNHelper.h" +#include "FWCore/MessageLogger/interface/MessageLogger.h" +#include "FWCore/Utilities/interface/FileInPath.h" +#include +#include +using namespace egammaTools; + +EgammaDNNHelper::EgammaDNNHelper(const DNNConfiguration& cfg, + const ModelSelector& modelSelector, + const std::vector& availableVars) + : cfg_(cfg), modelSelector_(modelSelector), nModels_(cfg_.modelsFiles.size()), graphDefs_(cfg_.modelsFiles.size()) { + initTensorFlowGraphs(); + initScalerFiles(availableVars); +} + +void EgammaDNNHelper::initTensorFlowGraphs() { + // load the graph definition + LogDebug("EgammaDNNHelper") << "Loading " << nModels_ << " graphs"; + size_t i = 0; + for (const auto& model_file : cfg_.modelsFiles) { + graphDefs_[i] = + std::unique_ptr(tensorflow::loadGraphDef(edm::FileInPath(model_file).fullPath())); + i++; + } +} + +std::vector EgammaDNNHelper::getSessions() const { + std::vector sessions; + LogDebug("EgammaDNNHelper") << "Starting " << nModels_ << " TF sessions"; + for (const auto& graphDef : graphDefs_) { + sessions.push_back(tensorflow::createSession(graphDef.get())); + } + LogDebug("EgammaDNNHelper") << "TF sessions started"; + return sessions; +} + +void EgammaDNNHelper::initScalerFiles(const std::vector& availableVars) { + for (const auto& scaler_file : cfg_.scalersFiles) { + // Parse scaler configuration + std::vector features; + std::ifstream inputfile_scaler{edm::FileInPath(scaler_file).fullPath()}; + int ninputs = 0; + if (inputfile_scaler.fail()) { + throw cms::Exception("MissingFile") << "Scaler file for Electron PFid DNN not found"; + } else { + // Now read mean, scale factors for each variable + float par1, par2; + std::string varName, type_str; + uint type; + while (inputfile_scaler >> varName >> type_str >> par1 >> par2) { + if (type_str == "stdscale") + type = 1; + else if (type_str == "minmax") + type = 2; + else if (type_str == "custom1") // 2*((X_train - minValues)/(MaxMinusMin)) -1.0 + type = 3; + else + type = 0; + features.push_back(ScalerConfiguration{.varName = varName, .type = type, .par1 = par1, .par2 = par2}); + // Protection for mismatch between requested variables and the available ones + auto match = std::find(availableVars.begin(), availableVars.end(), varName); + if (match == std::end(availableVars)) { + throw cms::Exception("MissingVariable") + << "Requested variable (" << varName << ") not available between DNN inputs"; + } + ninputs += 1; + } + } + inputfile_scaler.close(); + featuresMap_.push_back(features); + nInputs_.push_back(ninputs); + } +} + +std::pair> EgammaDNNHelper::getScaledInputs( + const std::map& variables) const { + // Call the modelSelector function passing the variables map to return + // the modelIndex to be used for the current candidate + const auto modelIndex = modelSelector_(variables); + std::vector inputs; + // Loop on the list of requested variables and scaling values for the specific modelIndex + // Different type of scaling are available: 0=no scaling, 1=standard scaler, 2=minmax + for (auto& [varName, type, par1, par2] : featuresMap_[modelIndex]) { + if (type == 1) // Standard scaling + inputs.push_back((variables.at(varName) - par1) / par2); + else if (type == 2) // MinMax + inputs.push_back((variables.at(varName) - par1) / (par2 - par1)); + else if (type == 3) //2*((X_train - minValues)/(MaxMinusMin)) -1.0 + inputs.push_back(2 * (variables.at(varName) - par1) / (par2 - par1) - 1.); + else { + inputs.push_back(variables.at(varName)); // Do nothing on the variable + } + //Protection for mismatch between requested variables and the available ones + // have been added when the scaler config are loaded --> here we know that the variables are available + } + return std::make_pair(modelIndex, inputs); +} + +std::vector> EgammaDNNHelper::evaluate(const std::vector>& candidates, + const std::vector& sessions) const { + /* + Evaluate the PFID DNN for all the electrons/photons. + nModels_ are defined depending on modelIndex --> we need to build N input tensors to evaluate + the DNNs with batching. + + 1) Get all the variable for each candidate vector> + 2) Scale the input and select the variables for each model + 2) Prepare the input tensors for the models + 3) Run the models and get the output for each candidate + 4) Sort the output by candidate index + 5) Return the DNN outputs + + */ + size_t nCandidates = candidates.size(); + std::vector> indexMap(nModels_); // for each model; the list of candidate index is saved + std::vector> inputsVectors(nCandidates); + std::vector counts(nModels_); + + LogDebug("EgammaDNNHelper") << "Working on " << nCandidates << " candidates"; + + int icand = 0; + for (auto& candidate : candidates) { + LogDebug("EgammaDNNHelper") << "Working on candidate: " << icand; + const auto& [model_index, inputs] = getScaledInputs(candidate); + counts[model_index] += 1; + indexMap[model_index].push_back(icand); + inputsVectors[icand] = inputs; + icand++; + } + + // Prepare one input tensors for each model + std::vector input_tensors(nModels_); + // Pointers for filling efficiently the input tensors + std::vector input_tensors_pointer(nModels_); + for (size_t i = 0; i < nModels_; i++) { + LogDebug("EgammaDNNHelper") << "Initializing TF input " << i << " with rows:" << counts[i] + << " and cols:" << nInputs_[i]; + input_tensors[i] = tensorflow::Tensor{tensorflow::DT_FLOAT, {counts[i], nInputs_[i]}}; + input_tensors_pointer[i] = input_tensors[i].flat().data(); + } + + // Filling the input tensors + for (size_t m = 0; m < nModels_; m++) { + LogDebug("EgammaDNNHelper") << "Loading TF input tensor for model: " << m; + float* T = input_tensors_pointer[m]; + for (size_t cand_index : indexMap[m]) { + for (size_t k = 0; k < nInputs_[m]; k++, T++) { //Note the input tensor pointer incremented + *T = inputsVectors[cand_index][k]; + } + } + } + + // Define the output and run + // Define the output and run + std::vector>> outputs; + // Run all the models + for (size_t m = 0; m < nModels_; m++) { + if (counts[m] == 0) + continue; //Skip model witout inputs + std::vector output; + LogDebug("EgammaDNNHelper") << "Run model: " << m << " with " << counts[m] << " electrons"; + tensorflow::run(sessions[m], {{cfg_.inputTensorName, input_tensors[m]}}, {cfg_.outputTensorName}, &output); + // Get the output and save the ElectronDNNEstimator::outputDim numbers along with the ele index + const auto& r = output[0].tensor(); + // Iterate on the list of elements in the batch --> many electrons + for (uint b = 0; b < counts[m]; b++) { + std::vector result(cfg_.outputDim); + for (size_t k = 0; k < cfg_.outputDim; k++) + result[k] = r(b, k); + // Get the original index of the electorn in the original order + const auto cand_index = indexMap[m][b]; + outputs.push_back(std::make_pair(cand_index, result)); + } + } + // Now we have just to re-order the outputs + std::sort(outputs.begin(), outputs.end()); + std::vector> final_outputs(outputs.size()); + std::transform(outputs.begin(), outputs.end(), final_outputs.begin(), [](auto a) { return a.second; }); + + return final_outputs; +} diff --git a/RecoEgamma/ElectronIdentification/BuildFile.xml b/RecoEgamma/ElectronIdentification/BuildFile.xml index 42e38970d1813..b32916ed44533 100644 --- a/RecoEgamma/ElectronIdentification/BuildFile.xml +++ b/RecoEgamma/ElectronIdentification/BuildFile.xml @@ -10,6 +10,7 @@ + diff --git a/RecoEgamma/ElectronIdentification/interface/ElectronDNNEstimator.h b/RecoEgamma/ElectronIdentification/interface/ElectronDNNEstimator.h new file mode 100644 index 0000000000000..aabbfe954efbc --- /dev/null +++ b/RecoEgamma/ElectronIdentification/interface/ElectronDNNEstimator.h @@ -0,0 +1,41 @@ +#ifndef RecoEgamma_ElectronIdentification_ElectronDNNEstimator_h +#define RecoEgamma_ElectronIdentification_ElectronDNNEstimator_h + +#include "RecoEgamma/EgammaTools/interface/EgammaDNNHelper.h" +#include "PhysicsTools/TensorFlow/interface/TensorFlow.h" +#include "DataFormats/EgammaCandidates/interface/GsfElectron.h" + +#include +#include +#include + +class ElectronDNNEstimator { +public: + ElectronDNNEstimator(const egammaTools::DNNConfiguration&, const bool useEBModelInGap); + + std::vector getSessions() const; + ; + + // Function returning a map with all the possible variables and their name + std::map getInputsVars(const reco::GsfElectron& ele) const; + + // Evaluate the DNN on all the electrons with the correct model + std::vector> evaluate(const reco::GsfElectronCollection& ele, + const std::vector& sessions) const; + + // List of input variables names used to check the variables request as + // inputs in a dynamic way from configuration file. + // If an input variables is not found at construction time an expection is thrown. + static const std::vector dnnAvaibleInputs; + + static constexpr float ptThreshold = 10.; + static constexpr float ecalBarrelMaxEtaWithGap = 1.566; + static constexpr float ecalBarrelMaxEtaNoGap = 1.485; + +private: + const egammaTools::EgammaDNNHelper dnnHelper_; + + const bool useEBModelInGap_; +}; + +#endif diff --git a/RecoEgamma/ElectronIdentification/src/ElectronDNNEstimator.cc b/RecoEgamma/ElectronIdentification/src/ElectronDNNEstimator.cc new file mode 100644 index 0000000000000..a9efb07025ac2 --- /dev/null +++ b/RecoEgamma/ElectronIdentification/src/ElectronDNNEstimator.cc @@ -0,0 +1,110 @@ +#include "RecoEgamma/ElectronIdentification/interface/ElectronDNNEstimator.h" +#include "FWCore/MessageLogger/interface/MessageLogger.h" +#include "FWCore/Utilities/interface/FileInPath.h" + +#include +#include +#include + +using namespace std::placeholders; + +inline uint electronModelSelector(const std::map& vars, float ptThr, float etaThr) { + /* + Selection of the model to be applied on the electron based on pt/eta cuts or whatever selection + */ + const auto pt = vars.at("pt"); + const auto absEta = std::abs(vars.at("eta")); + if (pt < ptThr) + return 0; + else { + if (absEta <= etaThr) { + return 1; + } else { + return 2; + } + } +} + +ElectronDNNEstimator::ElectronDNNEstimator(const egammaTools::DNNConfiguration& cfg, const bool useEBModelInGap) + : dnnHelper_(cfg, + std::bind(electronModelSelector, + _1, + ElectronDNNEstimator::ptThreshold, + (useEBModelInGap) ? ElectronDNNEstimator::ecalBarrelMaxEtaWithGap + : ElectronDNNEstimator::ecalBarrelMaxEtaNoGap), + ElectronDNNEstimator::dnnAvaibleInputs), + useEBModelInGap_(useEBModelInGap) {} + +std::vector ElectronDNNEstimator::getSessions() const { return dnnHelper_.getSessions(); }; + +const std::vector ElectronDNNEstimator::dnnAvaibleInputs = {{"pt", + "eta", + "fbrem", + "abs(deltaEtaSuperClusterTrackAtVtx)", + "abs(deltaPhiSuperClusterTrackAtVtx)", + "full5x5_sigmaIetaIeta", + "full5x5_hcalOverEcal", + "eSuperClusterOverP", + "full5x5_e1x5", + "eEleClusterOverPout", + "closestCtfTrackNormChi2", + "closestCtfTrackNLayers", + "gsfTrack.missing_inner_hits", + "dr03TkSumPt", + "dr03EcalRecHitSumEt", + "dr03HcalTowerSumEt", + "gsfTrack.normalizedChi2", + "superCluster.eta", + "ecalPFClusterIso", + "hcalPFClusterIso", + "numberOfBrems", + "abs(deltaEtaSeedClusterTrackAtCalo)", + "hadronicOverEm", + "full5x5_e2x5Max", + "full5x5_e5x5"}}; + +std::map ElectronDNNEstimator::getInputsVars(const reco::GsfElectron& ele) const { + // Prepare a map with all the defined variables + std::map variables; + reco::TrackRef myTrackRef = ele.closestCtfTrackRef(); + bool validKF = (myTrackRef.isNonnull() && myTrackRef.isAvailable()); + variables["pt"] = ele.pt(); + variables["eta"] = ele.eta(); + variables["fbrem"] = ele.fbrem(); + variables["abs(deltaEtaSuperClusterTrackAtVtx)"] = std::abs(ele.deltaEtaSuperClusterTrackAtVtx()); + variables["abs(deltaPhiSuperClusterTrackAtVtx)"] = std::abs(ele.deltaPhiSuperClusterTrackAtVtx()); + variables["full5x5_sigmaIetaIeta"] = ele.full5x5_sigmaIetaIeta(); + variables["full5x5_hcalOverEcal"] = ele.full5x5_hcalOverEcalValid() ? ele.full5x5_hcalOverEcal() : 0; + variables["eSuperClusterOverP"] = ele.eSuperClusterOverP(); + variables["full5x5_e1x5"] = ele.full5x5_e1x5(); + variables["eEleClusterOverPout"] = ele.eEleClusterOverPout(); + variables["closestCtfTrackNormChi2"] = ele.closestCtfTrackNormChi2(); + variables["closestCtfTrackNLayers"] = ele.closestCtfTrackNLayers(); + variables["gsfTrack.missing_inner_hits"] = + (validKF) ? myTrackRef->hitPattern().numberOfLostHits(reco::HitPattern::MISSING_INNER_HITS) : -1.; + variables["dr03TkSumPt"] = ele.dr03TkSumPt(); + variables["dr03EcalRecHitSumEt"] = ele.dr03EcalRecHitSumEt(); + variables["dr03HcalTowerSumEt"] = ele.dr03HcalTowerSumEt(); + variables["gsfTrack.normalizedChi2"] = (validKF) ? myTrackRef->normalizedChi2() : 0; + variables["superCluster.eta"] = ele.superCluster()->eta(); + variables["ecalPFClusterIso"] = ele.ecalPFClusterIso(); + variables["hcalPFClusterIso"] = ele.hcalPFClusterIso(); + variables["numberOfBrems"] = ele.numberOfBrems(); + variables["abs(deltaEtaSeedClusterTrackAtCalo)"] = std::abs(ele.deltaEtaSeedClusterTrackAtCalo()); + variables["hadronicOverEm"] = ele.hcalOverEcalValid() ? ele.hadronicOverEm() : 0; + variables["full5x5_e2x5Max"] = ele.full5x5_e2x5Max(); + variables["full5x5_e5x5"] = ele.full5x5_e5x5(); + // Define more variables here and use them directly in the model config! + return variables; +} + +std::vector> ElectronDNNEstimator::evaluate( + const reco::GsfElectronCollection& electrons, const std::vector& sessions) const { + // Collect the map of variables for each candidate and call the dnnHelper + // Scaling, model selection and running is performed in the helper + std::vector> inputs; + for (const auto& ele : electrons) { + inputs.push_back(getInputsVars(ele)); + } + return dnnHelper_.evaluate(inputs, sessions); +} diff --git a/RecoEgamma/PhotonIdentification/BuildFile.xml b/RecoEgamma/PhotonIdentification/BuildFile.xml index df2125cded3cc..02e53c5fd2adf 100644 --- a/RecoEgamma/PhotonIdentification/BuildFile.xml +++ b/RecoEgamma/PhotonIdentification/BuildFile.xml @@ -20,8 +20,10 @@ + + diff --git a/RecoEgamma/PhotonIdentification/interface/PhotonDNNEstimator.h b/RecoEgamma/PhotonIdentification/interface/PhotonDNNEstimator.h new file mode 100644 index 0000000000000..d07aff1c932e8 --- /dev/null +++ b/RecoEgamma/PhotonIdentification/interface/PhotonDNNEstimator.h @@ -0,0 +1,41 @@ +#ifndef RecoEgamma_PhotonIdentification_PhotonDNNEstimator_h +#define RecoEgamma_PhotonIdentification_PhotonDNNEstimator_h + +#include "RecoEgamma/EgammaTools/interface/EgammaDNNHelper.h" +#include "PhysicsTools/TensorFlow/interface/TensorFlow.h" +#include "DataFormats/EgammaCandidates/interface/Photon.h" +#include "DataFormats/EgammaCandidates/interface/PhotonFwd.h" + +#include +#include +#include + +class PhotonDNNEstimator { +public: + PhotonDNNEstimator(const egammaTools::DNNConfiguration&, const bool useEBModelInGap); + + std::vector getSessions() const; + ; + + // Function returning a map with all the possible variables and their name + std::map getInputsVars(const reco::Photon& ele) const; + + // Evaluate the DNN on all the electrons with the correct model + std::vector> evaluate(const reco::PhotonCollection& ele, + const std::vector& sessions) const; + + // List of input variables names used to check the variables request as + // inputs in a dynamic way from configuration file. + // If an input variables is not found at construction time an expection is thrown. + static const std::vector dnnAvaibleInputs; + + static constexpr float ecalBarrelMaxEtaWithGap = 1.566; + static constexpr float ecalBarrelMaxEtaNoGap = 1.485; + +private: + const egammaTools::EgammaDNNHelper dnnHelper_; + + const bool useEBModelInGap_; +}; + +#endif diff --git a/RecoEgamma/PhotonIdentification/src/PhotonDNNEstimator.cc b/RecoEgamma/PhotonIdentification/src/PhotonDNNEstimator.cc new file mode 100644 index 0000000000000..0207282ad822a --- /dev/null +++ b/RecoEgamma/PhotonIdentification/src/PhotonDNNEstimator.cc @@ -0,0 +1,78 @@ +#include "RecoEgamma/PhotonIdentification/interface/PhotonDNNEstimator.h" +#include "FWCore/MessageLogger/interface/MessageLogger.h" +#include "FWCore/Utilities/interface/FileInPath.h" + +#include +#include +#include + +using namespace std::placeholders; + +inline uint photonModelSelector(const std::map& vars, float etaThr) { + /* + Selection of the model to be applied on the photon based on eta limit + */ + const auto absEta = std::abs(vars.at("eta")); + if (absEta <= etaThr) { + return 0; + } else { + return 1; + } +} + +PhotonDNNEstimator::PhotonDNNEstimator(const egammaTools::DNNConfiguration& cfg, const bool useEBModelInGap) + : dnnHelper_(cfg, + std::bind(photonModelSelector, + _1, + (useEBModelInGap) ? PhotonDNNEstimator::ecalBarrelMaxEtaWithGap + : PhotonDNNEstimator::ecalBarrelMaxEtaNoGap), + PhotonDNNEstimator::dnnAvaibleInputs), + useEBModelInGap_(useEBModelInGap) {} + +std::vector PhotonDNNEstimator::getSessions() const { return dnnHelper_.getSessions(); }; + +const std::vector PhotonDNNEstimator::dnnAvaibleInputs = {{"pt", + "eta", + "hadTowOverEm", + "TrkSumPtHollow", + "EcalRecHit", + "SigmaIetaIeta", + "SigmaIEtaIEtaFull5x5", + "SigmaIEtaIPhiFull5x5", + "EcalPFClusterIso", + "HcalPFClusterIso", + "HasPixelSeed", + "R9Full5x5", + "hcalTower"}}; + +std::map PhotonDNNEstimator::getInputsVars(const reco::Photon& photon) const { + // Prepare a map with all the defined variables + std::map variables; + variables["pt"] = photon.pt(); + variables["eta"] = photon.eta(); + variables["hadTowOverEm"] = photon.hadTowOverEmValid() ? photon.hadTowOverEm() : 0; + variables["TrkSumPtHollow"] = photon.trkSumPtHollowConeDR03(); + variables["EcalRecHit"] = photon.ecalRecHitSumEtConeDR03(); + variables["SigmaIetaIeta"] = photon.sigmaIetaIeta(); + variables["SigmaIEtaIEtaFull5x5"] = photon.full5x5_sigmaIetaIeta(); + variables["SigmaIEtaIPhiFull5x5"] = photon.full5x5_showerShapeVariables().sigmaIetaIphi; + variables["EcalPFClusterIso"] = photon.ecalPFClusterIso(); + variables["HcalPFClusterIso"] = photon.hcalPFClusterIso(); + variables["HasPixelSeed"] = (Int_t)photon.hasPixelSeed(); + variables["R9Full5x5"] = photon.full5x5_r9(); + variables["hcalTower"] = photon.hcalTowerSumEtConeDR03(); + variables["R9Full5x5"] = photon.full5x5_r9(); + // Define more variables here and use them directly in the model config! + return variables; +} + +std::vector> PhotonDNNEstimator::evaluate(const reco::PhotonCollection& photons, + const std::vector& sessions) const { + // Collect the map of variables for each candidate and call the dnnHelper + // Scaling, model selection and running is performed in the helper + std::vector> inputs; + for (const auto& photon : photons) { + inputs.push_back(getInputsVars(photon)); + } + return dnnHelper_.evaluate(inputs, sessions); +} diff --git a/RecoParticleFlow/PFProducer/interface/PFEGammaFilters.h b/RecoParticleFlow/PFProducer/interface/PFEGammaFilters.h index 234c90e52ddc5..27861d8c1fc3e 100644 --- a/RecoParticleFlow/PFProducer/interface/PFEGammaFilters.h +++ b/RecoParticleFlow/PFProducer/interface/PFEGammaFilters.h @@ -44,12 +44,23 @@ class PFEGammaFilters { float pho_sumPtTrackIsoSlope_; // Electron selections + const bool useElePFidDNN_; + const bool usePhotonPFidDNN_; + const bool useEBModelInGap_; const float ele_iso_pt_; const float ele_iso_mva_eb_; const float ele_iso_mva_ee_; const float ele_iso_combIso_eb_; const float ele_iso_combIso_ee_; const float ele_noniso_mva_; + // Threshold for DNN ele pfid + float ele_dnnLowPtThr_; + float ele_dnnHighPtBarrelThr_; + float ele_dnnHighPtEndcapThr_; + // Threshold for DNN photon pfid + float photon_dnnBarrelThr_; + float photon_dnnEndcapThr_; + const int ele_missinghits_; const float ele_ecalDrivenHademPreselCut_; const float ele_maxElePtForOnlyMVAPresel_; diff --git a/RecoParticleFlow/PFProducer/python/particleFlow_cff.py b/RecoParticleFlow/PFProducer/python/particleFlow_cff.py index 24f0272ccd63b..944b01ac41989 100644 --- a/RecoParticleFlow/PFProducer/python/particleFlow_cff.py +++ b/RecoParticleFlow/PFProducer/python/particleFlow_cff.py @@ -14,6 +14,17 @@ ## In 12_2, we expect to have EEE's ID/regression, then this switch can flip to True particleFlowTmp.PFEGammaFiltersParameters.allowEEEinPF = cms.bool(False) +# Thresholds for e/gamma PFID DNN +particleFlowTmp.PFEGammaFiltersParameters.electronDnnThresholds = cms.PSet( + electronDnnHighPtBarrelThr = cms.double(0.122), + electronDnnHighPtEndcapThr = cms.double(0.072), + electronDnnLowPtThr = cms.double(0.081) + ) +particleFlowTmp.PFEGammaFiltersParameters.photonDnnThresholds = cms.PSet( + photonDnnBarrelThr = cms.double(0.70), + photonDnnEndcapThr = cms.double(0.79) +) + from Configuration.Eras.Modifier_pf_badHcalMitigationOff_cff import pf_badHcalMitigationOff pf_badHcalMitigationOff.toModify(particleFlowTmp.PFEGammaFiltersParameters, electron_protectionsForBadHcal = dict(enableProtections = False), @@ -21,3 +32,10 @@ from Configuration.ProcessModifiers.egamma_lowPt_exclusive_cff import egamma_lowPt_exclusive egamma_lowPt_exclusive.toModify(particleFlowTmp.PFEGammaFiltersParameters,photon_MinEt = 1.) + +# Activate Egamma PFID with DNN for Run3 +from Configuration.Eras.Modifier_run3_common_cff import run3_common +run3_common.toModify(particleFlowTmp.PFEGammaFiltersParameters, + useElePFidDnn = True, + usePhotonPFidDnn = True +) diff --git a/RecoParticleFlow/PFProducer/src/PFAlgo.cc b/RecoParticleFlow/PFProducer/src/PFAlgo.cc index 5fb624b0a046b..94e887a788577 100644 --- a/RecoParticleFlow/PFProducer/src/PFAlgo.cc +++ b/RecoParticleFlow/PFProducer/src/PFAlgo.cc @@ -294,6 +294,12 @@ void PFAlgo::egammaFilters(const reco::PFBlockRef& blockref, myPFElectron.set_mva_e_pi(gedEleRef->mva_e_pi()); myPFElectron.set_mva_Isolated(gedEleRef->mva_Isolated()); + myPFElectron.set_dnn_e_sigIsolated(gedEleRef->dnn_signal_Isolated()); + myPFElectron.set_dnn_e_sigNonIsolated(gedEleRef->dnn_signal_nonIsolated()); + myPFElectron.set_dnn_e_bkgNonIsolated(gedEleRef->dnn_bkg_nonIsolated()); + myPFElectron.set_dnn_e_bkgTau(gedEleRef->dnn_bkg_Tau()); + myPFElectron.set_dnn_e_bkgPhoton(gedEleRef->dnn_bkg_Photon()); + LogTrace("PFAlgo|egammaFilters") << " PFAlgo: found an electron with NEW EGamma code "; LogTrace("PFAlgo|egammaFilters") << " myPFElectron: pt " << myPFElectron.pt() << " eta,phi " << myPFElectron.eta() << ", " << myPFElectron.phi() << " mva " @@ -339,6 +345,8 @@ void PFAlgo::egammaFilters(const reco::PFBlockRef& blockref, ::math::XYZPoint v(primaryVertex_.x(), primaryVertex_.y(), primaryVertex_.z()); myPFPhoton.setVertex(v); myPFPhoton.setP4(gedPhoRef->p4()); + // DNN pfid + myPFPhoton.set_dnn_gamma(gedPhoRef->pfDNN()); LogTrace("PFAlgo|egammaFilters") << " PFAlgo: found a photon with NEW EGamma code "; LogTrace("PFAlgo|egammaFilters") << " myPFPhoton: pt " << myPFPhoton.pt() << " eta,phi " << myPFPhoton.eta() << ", " << myPFPhoton.phi() << " charge " << myPFPhoton.charge(); diff --git a/RecoParticleFlow/PFProducer/src/PFEGammaFilters.cc b/RecoParticleFlow/PFProducer/src/PFEGammaFilters.cc index e5bcf510973e3..ad770b07eb62e 100644 --- a/RecoParticleFlow/PFProducer/src/PFEGammaFilters.cc +++ b/RecoParticleFlow/PFProducer/src/PFEGammaFilters.cc @@ -12,6 +12,10 @@ using namespace reco; namespace { + // Constants defining the ECAL barrel limit + constexpr float ecalBarrelMaxEtaWithGap = 1.566; + constexpr float ecalBarrelMaxEtaNoGap = 1.485; + void readEBEEParams_(const edm::ParameterSet& pset, const std::string& name, std::array& out) { const auto& vals = pset.getParameter>(name); if (vals.size() != 2) @@ -28,6 +32,9 @@ PFEGammaFilters::PFEGammaFilters(const edm::ParameterSet& cfg) ph_loose_hoe_(cfg.getParameter("photon_HoE")), ph_sietaieta_eb_(cfg.getParameter("photon_SigmaiEtaiEta_barrel")), ph_sietaieta_ee_(cfg.getParameter("photon_SigmaiEtaiEta_endcap")), + useElePFidDNN_(cfg.getParameter("useElePFidDnn")), + usePhotonPFidDNN_(cfg.getParameter("usePhotonPFidDnn")), + useEBModelInGap_(cfg.getParameter("useEBModelInGap")), ele_iso_pt_(cfg.getParameter("electron_iso_pt")), ele_iso_mva_eb_(cfg.getParameter("electron_iso_mva_barrel")), ele_iso_mva_ee_(cfg.getParameter("electron_iso_mva_endcap")), @@ -42,6 +49,8 @@ PFEGammaFilters::PFEGammaFilters(const edm::ParameterSet& cfg) auto const& eleProtectionsForJetMET = cfg.getParameter("electron_protectionsForJetMET"); auto const& phoProtectionsForBadHcal = cfg.getParameter("photon_protectionsForBadHcal"); auto const& phoProtectionsForJetMET = cfg.getParameter("photon_protectionsForJetMET"); + auto const& eleDNNIdThresholds = cfg.getParameter("electronDnnThresholds"); + auto const& photonDNNIdThresholds = cfg.getParameter("photonDnnThresholds"); pho_sumPtTrackIso_ = phoProtectionsForJetMET.getParameter("sumPtTrackIso"); pho_sumPtTrackIsoSlope_ = phoProtectionsForJetMET.getParameter("sumPtTrackIsoSlope"); @@ -60,6 +69,12 @@ PFEGammaFilters::PFEGammaFilters(const edm::ParameterSet& cfg) ele_maxEeleOverPout_ = eleProtectionsForJetMET.getParameter("maxEeleOverPout"); ele_maxDPhiIN_ = eleProtectionsForJetMET.getParameter("maxDPhiIN"); + ele_dnnLowPtThr_ = eleDNNIdThresholds.getParameter("electronDnnLowPtThr"); + ele_dnnHighPtBarrelThr_ = eleDNNIdThresholds.getParameter("electronDnnHighPtBarrelThr"); + ele_dnnHighPtEndcapThr_ = eleDNNIdThresholds.getParameter("electronDnnHighPtEndcapThr"); + photon_dnnBarrelThr_ = photonDNNIdThresholds.getParameter("photonDnnBarrelThr"); + photon_dnnEndcapThr_ = photonDNNIdThresholds.getParameter("photonDnnEndcapThr"); + readEBEEParams_(eleProtectionsForBadHcal, "full5x5_sigmaIetaIeta", badHcal_full5x5_sigmaIetaIeta_); readEBEEParams_(eleProtectionsForBadHcal, "eInvPInv", badHcal_eInvPInv_); readEBEEParams_(eleProtectionsForBadHcal, "dEta", badHcal_dEta_); @@ -92,26 +107,40 @@ bool PFEGammaFilters::passPhotonSelection(const reco::Photon& photon) const { : badHcal_phoTrkSolidConeIso_offs_ + badHcal_phoTrkSolidConeIso_slope_ * photon.pt()) << ")" << std::endl; - if (photon.hadTowOverEm() > ph_loose_hoe_) - return false; - //Isolation variables in 0.3 cone combined - if (photon.trkSumPtHollowConeDR03() + photon.ecalRecHitSumEtConeDR03() + photon.hcalTowerSumEtConeDR03() > - ph_combIso_) - return false; - - //patch for bad hcal - if (!validHoverE && badHcal_phoEnable_ && - photon.trkSumPtSolidConeDR03() > - badHcal_phoTrkSolidConeIso_offs_ + badHcal_phoTrkSolidConeIso_slope_ * photon.pt()) { - return false; - } - - if (photon.isEB()) { - if (photon.sigmaIetaIeta() > ph_sietaieta_eb_) - return false; + if (usePhotonPFidDNN_) { + // Run3 DNN based PFID + const auto dnn = photon.pfDNN(); + const auto photEta = std::abs(photon.eta()); + const auto etaThreshold = (useEBModelInGap_) ? ecalBarrelMaxEtaWithGap : ecalBarrelMaxEtaNoGap; + // using the Barrel model for photons in the EB-EE gap + if (photEta <= etaThreshold) { + return dnn > photon_dnnBarrelThr_; + } else if (photEta > etaThreshold) { + return dnn > photon_dnnEndcapThr_; + } } else { - if (photon.sigmaIetaIeta() > ph_sietaieta_ee_) + // Run2 cut based PFID + if (photon.hadTowOverEm() > ph_loose_hoe_) + return false; + //Isolation variables in 0.3 cone combined + if (photon.trkSumPtHollowConeDR03() + photon.ecalRecHitSumEtConeDR03() + photon.hcalTowerSumEtConeDR03() > + ph_combIso_) + return false; + + //patch for bad hcal + if (!validHoverE && badHcal_phoEnable_ && + photon.trkSumPtSolidConeDR03() > + badHcal_phoTrkSolidConeIso_offs_ + badHcal_phoTrkSolidConeIso_slope_ * photon.pt()) { return false; + } + + if (photon.isEB()) { + if (photon.sigmaIetaIeta() > ph_sietaieta_eb_) + return false; + } else { + if (photon.sigmaIetaIeta() > ph_sietaieta_ee_) + return false; + } } return true; @@ -136,31 +165,50 @@ bool PFEGammaFilters::passElectronSelection(const reco::GsfElectron& electron, bool passEleSelection = false; // Electron ET - float electronPt = electron.pt(); - - if (electronPt > ele_iso_pt_) { - double isoDr03 = electron.dr03TkSumPt() + electron.dr03EcalRecHitSumEt() + electron.dr03HcalTowerSumEt(); - double eleEta = fabs(electron.eta()); - if (eleEta <= 1.485 && isoDr03 < ele_iso_combIso_eb_) { - if (electron.mva_Isolated() > ele_iso_mva_eb_) - passEleSelection = true; - } else if (eleEta > 1.485 && isoDr03 < ele_iso_combIso_ee_) { - if (electron.mva_Isolated() > ele_iso_mva_ee_) - passEleSelection = true; + const auto electronPt = electron.pt(); + const auto eleEta = std::abs(electron.eta()); + + if (useElePFidDNN_) { // Use DNN for ele pfID >=CMSSW12_1 + const auto dnn_sig = electron.dnn_signal_Isolated() + electron.dnn_signal_nonIsolated(); + const auto etaThreshold = (useEBModelInGap_) ? ecalBarrelMaxEtaWithGap : ecalBarrelMaxEtaNoGap; + if (electronPt > ele_iso_pt_) { + // using the Barrel model for electron in the EB-EE gap + if (eleEta <= etaThreshold) { + passEleSelection = dnn_sig > ele_dnnHighPtBarrelThr_; + } else if (eleEta > etaThreshold) { + passEleSelection = dnn_sig > ele_dnnHighPtEndcapThr_; + } + } else { // pt < ele_iso_pt_ + passEleSelection = dnn_sig > ele_dnnLowPtThr_; + } + // TODO: For the moment do not evaluate further conditions on isolation and HCAL cleaning.. + // To be understood if they are needed + + } else { // Use legacy MVA for ele pfID < CMSSW_12_1 + if (electronPt > ele_iso_pt_) { + double isoDr03 = electron.dr03TkSumPt() + electron.dr03EcalRecHitSumEt() + electron.dr03HcalTowerSumEt(); + if (eleEta <= ecalBarrelMaxEtaNoGap && isoDr03 < ele_iso_combIso_eb_) { + if (electron.mva_Isolated() > ele_iso_mva_eb_) + passEleSelection = true; + } else if (eleEta > ecalBarrelMaxEtaNoGap && isoDr03 < ele_iso_combIso_ee_) { + if (electron.mva_Isolated() > ele_iso_mva_ee_) + passEleSelection = true; + } } - } - // cout << " My OLD MVA " << pfcand.mva_e_pi() << " MyNEW MVA " << electron.mva() << endl; - if (electron.mva_e_pi() > ele_noniso_mva_) { - if (validHoverE || !badHcal_eleEnable_) { - passEleSelection = true; - } else { - bool EE = (std::abs(electron.eta()) > 1.485); // for prefer consistency with above than with E/gamma for now - if ((electron.full5x5_sigmaIetaIeta() < badHcal_full5x5_sigmaIetaIeta_[EE]) && - (std::abs(1.0 - electron.eSuperClusterOverP()) / electron.ecalEnergy() < badHcal_eInvPInv_[EE]) && - (std::abs(electron.deltaEtaSeedClusterTrackAtVtx()) < badHcal_dEta_[EE]) && // looser in case of misalignment - (std::abs(electron.deltaPhiSuperClusterTrackAtVtx()) < badHcal_dPhi_[EE])) { + if (electron.mva_e_pi() > ele_noniso_mva_) { + if (validHoverE || !badHcal_eleEnable_) { passEleSelection = true; + } else { + bool EE = (std::abs(electron.eta()) > + ecalBarrelMaxEtaNoGap); // for prefer consistency with above than with E/gamma for now + if ((electron.full5x5_sigmaIetaIeta() < badHcal_full5x5_sigmaIetaIeta_[EE]) && + (std::abs(1.0 - electron.eSuperClusterOverP()) / electron.ecalEnergy() < badHcal_eInvPInv_[EE]) && + (std::abs(electron.deltaEtaSeedClusterTrackAtVtx()) < + badHcal_dEta_[EE]) && // looser in case of misalignment + (std::abs(electron.deltaPhiSuperClusterTrackAtVtx()) < badHcal_dPhi_[EE])) { + passEleSelection = true; + } } } } @@ -293,10 +341,11 @@ bool PFEGammaFilters::isElectronSafeForJetMET(const reco::GsfElectron& electron, isSafeForJetMET = false; } // the electron is retained and the kf tracks are not locked - if ((fabs(1. - EtotPinMode) < ele_maxEcalEOverPRes_ && - (fabs(electron.eta()) < 1.0 || fabs(electron.eta()) > 2.0)) || - ((EtotPinMode < 1.1 && EtotPinMode > 0.6) && (fabs(electron.eta()) >= 1.0 && fabs(electron.eta()) <= 2.0))) { - if (fabs(1. - EGsfPoutMode) < ele_maxEeleOverPoutRes_ && (itrackHcalLinked == iextratrack)) { + if ((std::abs(1. - EtotPinMode) < ele_maxEcalEOverPRes_ && + (std::abs(electron.eta()) < 1.0 || std::abs(electron.eta()) > 2.0)) || + ((EtotPinMode < 1.1 && EtotPinMode > 0.6) && + (std::abs(electron.eta()) >= 1.0 && std::abs(electron.eta()) <= 2.0))) { + if (std::abs(1. - EGsfPoutMode) < ele_maxEeleOverPoutRes_ && (itrackHcalLinked == iextratrack)) { lockTracks = false; // lockExtraKf = false; if (debugSafeForJetMET) @@ -325,7 +374,7 @@ bool PFEGammaFilters::isElectronSafeForJetMET(const reco::GsfElectron& electron, } // For not-preselected Gsf Tracks ET > 50 GeV, apply dphi preselection - if (ETtotal > ele_maxE_ && fabs(dphi_normalsc) > ele_maxDPhiIN_) { + if (ETtotal > ele_maxE_ && std::abs(dphi_normalsc) > ele_maxDPhiIN_) { if (debugSafeForJetMET) cout << " *****This electron candidate is discarded Large ANGLE " << " ETtotal " << ETtotal << " EGsfPoutMode " << dphi_normalsc << endl; @@ -441,6 +490,24 @@ void PFEGammaFilters::fillPSetDescription(edm::ParameterSetDescription& iDesc) { iDesc.add("electron_ecalDrivenHademPreselCut", 0.15); iDesc.add("electron_maxElePtForOnlyMVAPresel", 50.0); iDesc.add("allowEEEinPF", false); + iDesc.add("useElePFidDnn", false); + { + edm::ParameterSetDescription psd; + psd.add("electronDnnLowPtThr", 0.5); + psd.add("electronDnnHighPtBarrelThr", 0.5); + psd.add("electronDnnHighPtEndcapThr", 0.5); + + iDesc.add("electronDnnThresholds", psd); + } + iDesc.add("usePhotonPFidDnn", false); + { + edm::ParameterSetDescription psd; + psd.add("photonDnnBarrelThr", 0.5); + psd.add("photonDnnEndcapThr", 0.5); + iDesc.add("photonDnnThresholds", psd); + } + // control if the EB DNN models should be used up to eta 1.485 or 1.566 + iDesc.add("useEBModelInGap", true); { edm::ParameterSetDescription psd; psd.add("maxNtracks", 3.0)->setComment("Max tracks pointing at Ele cluster");