File indexing completed on 2022-06-16 01:50:14
0001 #include "CommonTools/Utils/interface/StringToEnumValue.h"
0002 #include "DataFormats/Common/interface/Handle.h"
0003 #include "DataFormats/EcalRecHit/interface/EcalSeverityLevel.h"
0004 #include "DataFormats/EgammaCandidates/interface/GsfElectron.h"
0005 #include "DataFormats/EgammaReco/interface/ElectronSeed.h"
0006 #include "DataFormats/GsfTrackReco/interface/GsfTrack.h"
0007 #include "DataFormats/ParticleFlowCandidate/interface/PFCandidate.h"
0008 #include "DataFormats/ParticleFlowCandidate/interface/PFCandidateEGammaExtra.h"
0009 #include "DataFormats/ParticleFlowReco/interface/GsfPFRecTrack.h"
0010 #include "DataFormats/TrackCandidate/interface/TrackCandidateCollection.h"
0011 #include "FWCore/Common/interface/Provenance.h"
0012 #include "FWCore/Framework/interface/Event.h"
0013 #include "FWCore/Framework/interface/EventSetup.h"
0014 #include "FWCore/Framework/interface/stream/EDProducer.h"
0015 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0016 #include "FWCore/ParameterSet/interface/ConfigurationDescriptions.h"
0017 #include "FWCore/ParameterSet/interface/ParameterSet.h"
0018 #include "FWCore/ParameterSet/interface/ParameterSetDescription.h"
0019 #include "RecoEcal/EgammaCoreTools/interface/EcalClusterFunctionFactory.h"
0020 #include "RecoEgamma/EgammaIsolationAlgos/interface/EleTkIsolFromCands.h"
0021 #include "RecoEgamma/EgammaElectronAlgos/interface/EgAmbiguityTools.h"
0022 #include "RecoEgamma/EgammaElectronAlgos/interface/ElectronUtilities.h"
0023 #include "RecoEgamma/EgammaElectronAlgos/interface/GsfElectronAlgo.h"
0024 #include "RecoEcal/EgammaCoreTools/interface/EgammaLocalCovParamDefaults.h"
0025 #include "RecoEgamma/EgammaIsolationAlgos/interface/EcalPFClusterIsolation.h"
0026 #include "RecoEgamma/EgammaIsolationAlgos/interface/HcalPFClusterIsolation.h"
0027
0028 #include <array>
0029
0030 using namespace reco;
0031
0032 namespace {
0033
0034 void setMVAOutputs(reco::GsfElectronCollection& electrons,
0035 const GsfElectronAlgo::HeavyObjectCache* hoc,
0036 reco::VertexCollection const& vertices,
0037 bool dnnPFidEnabled,
0038 float extetaboundary,
0039 const std::vector<tensorflow::Session*>& tfSessions) {
0040 std::vector<GsfElectron::MvaOutput> mva_outputs(electrons.size());
0041 size_t iele = 0;
0042 for (auto& el : electrons) {
0043 GsfElectron::MvaOutput mvaOutput;
0044 mvaOutput.mva_e_pi = hoc->sElectronMVAEstimator->mva(el, vertices);
0045 mvaOutput.mva_Isolated = hoc->iElectronMVAEstimator->mva(el, vertices.size());
0046 if (dnnPFidEnabled) {
0047 mva_outputs[iele] = mvaOutput;
0048 } else {
0049 el.setMvaOutput(mvaOutput);
0050 }
0051 iele++;
0052 }
0053 if (dnnPFidEnabled) {
0054
0055 LogDebug("GsfElectronProducer") << "Getting DNN PFId for ele";
0056 const auto& dnn_ele_pfid = hoc->iElectronDNNEstimator->evaluate(electrons, tfSessions);
0057 int jele = 0;
0058 for (auto& el : electrons) {
0059 const auto& [iModel, values] = dnn_ele_pfid[jele];
0060
0061 auto& mvaOutput = mva_outputs[jele];
0062
0063 if (iModel <= 3) {
0064 assert(values.size() == 5);
0065 mvaOutput.dnn_e_sigIsolated = values[0];
0066 mvaOutput.dnn_e_sigNonIsolated = values[1];
0067 mvaOutput.dnn_e_bkgNonIsolated = values[2];
0068 mvaOutput.dnn_e_bkgTau = values[3];
0069 mvaOutput.dnn_e_bkgPhoton = values[4];
0070 } else if (iModel == 4) {
0071 assert(values.size() == 3);
0072 mvaOutput.dnn_e_sigIsolated = values[0];
0073 mvaOutput.dnn_e_sigNonIsolated = 0.0;
0074 mvaOutput.dnn_e_bkgNonIsolated = values[1];
0075 mvaOutput.dnn_e_bkgTau = 0.0;
0076 mvaOutput.dnn_e_bkgPhoton = values[2];
0077 }
0078 el.setMvaOutput(mvaOutput);
0079 jele++;
0080 }
0081 }
0082 }
0083
0084
0085
0086
0087 auto matchWithPFCandidates(std::vector<reco::PFCandidate> const& pfCandidates) {
0088 std::map<reco::GsfTrackRef, reco::GsfElectron::MvaInput> gsfMVAInputs{};
0089
0090
0091 for (auto const& pfCand : pfCandidates) {
0092 reco::GsfElectronRef myRef;
0093
0094 if (pfCand.gsfTrackRef().isNonnull()) {
0095 reco::GsfElectron::MvaInput input;
0096 input.earlyBrem = pfCand.egammaExtraRef()->mvaVariable(reco::PFCandidateEGammaExtra::MVA_FirstBrem);
0097 input.lateBrem = pfCand.egammaExtraRef()->mvaVariable(reco::PFCandidateEGammaExtra::MVA_LateBrem);
0098 input.deltaEta = pfCand.egammaExtraRef()->mvaVariable(reco::PFCandidateEGammaExtra::MVA_DeltaEtaTrackCluster);
0099 input.sigmaEtaEta = pfCand.egammaExtraRef()->sigmaEtaEta();
0100 input.hadEnergy = pfCand.egammaExtraRef()->hadEnergy();
0101 gsfMVAInputs[pfCand.gsfTrackRef()] = input;
0102 }
0103 }
0104 return gsfMVAInputs;
0105 }
0106
0107 void logElectrons(reco::GsfElectronCollection const& electrons, edm::Event const& event, const std::string& title) {
0108 LogTrace("GsfElectronAlgo") << "========== " << title << " ==========";
0109 LogTrace("GsfElectronAlgo") << "Event: " << event.id();
0110 LogTrace("GsfElectronAlgo") << "Number of electrons: " << electrons.size();
0111 for (auto const& ele : electrons) {
0112 LogTrace("GsfElectronAlgo") << "Electron with charge, pt, eta, phi: " << ele.charge() << " , " << ele.pt()
0113 << " , " << ele.eta() << " , " << ele.phi();
0114 }
0115 LogTrace("GsfElectronAlgo") << "=================================================";
0116 }
0117
0118 }
0119
0120 class GsfElectronProducer : public edm::stream::EDProducer<edm::GlobalCache<GsfElectronAlgo::HeavyObjectCache>> {
0121 public:
0122 static void fillDescriptions(edm::ConfigurationDescriptions&);
0123
0124 explicit GsfElectronProducer(const edm::ParameterSet&, const GsfElectronAlgo::HeavyObjectCache*);
0125
0126 static std::unique_ptr<GsfElectronAlgo::HeavyObjectCache> initializeGlobalCache(const edm::ParameterSet& conf) {
0127 return std::make_unique<GsfElectronAlgo::HeavyObjectCache>(conf);
0128 }
0129
0130 void endStream() override;
0131
0132 static void globalEndJob(GsfElectronAlgo::HeavyObjectCache const*){};
0133
0134
0135 void produce(edm::Event& event, const edm::EventSetup& setup) override;
0136
0137 private:
0138 std::unique_ptr<GsfElectronAlgo> algo_;
0139
0140
0141 GsfElectronAlgo::Tokens inputCfg_;
0142 GsfElectronAlgo::StrategyConfiguration strategyCfg_;
0143 const GsfElectronAlgo::CutsConfiguration cutsCfg_;
0144 ElectronHcalHelper::Configuration hcalCfg_, hcalCfgBc_;
0145
0146 bool hcalRun2EffDepth_;
0147
0148 bool isPreselected(reco::GsfElectron const& ele) const;
0149 void setAmbiguityData(reco::GsfElectronCollection& electrons,
0150 edm::Event const& event,
0151 bool ignoreNotPreselected = true) const;
0152
0153
0154 bool ecalSeedingParametersChecked_;
0155 void checkEcalSeedingParameters(edm::ParameterSet const&);
0156
0157 const edm::EDPutTokenT<reco::GsfElectronCollection> electronPutToken_;
0158 const edm::EDGetTokenT<reco::GsfPFRecTrackCollection> gsfPfRecTracksTag_;
0159 edm::EDGetTokenT<reco::PFCandidateCollection> egmPFCandidateCollection_;
0160
0161 const bool useGsfPfRecTracks_;
0162
0163 const bool resetMvaValuesUsingPFCandidates_;
0164
0165 bool dnnPFidEnabled_;
0166 float extetaboundary_;
0167
0168 std::vector<tensorflow::Session*> tfSessions_;
0169 };
0170
0171 void GsfElectronProducer::fillDescriptions(edm::ConfigurationDescriptions& descriptions) {
0172 edm::ParameterSetDescription desc;
0173
0174 desc.add<edm::InputTag>("gsfElectronCoresTag", {"ecalDrivenGsfElectronCores"});
0175 desc.add<edm::InputTag>("vtxTag", {"offlinePrimaryVertices"});
0176 desc.add<edm::InputTag>("conversionsTag", {"allConversions"});
0177 desc.add<edm::InputTag>("gsfPfRecTracksTag", {"pfTrackElec"});
0178 desc.add<edm::InputTag>("barrelRecHitCollectionTag", {"ecalRecHit", "EcalRecHitsEB"});
0179 desc.add<edm::InputTag>("endcapRecHitCollectionTag", {"ecalRecHit", "EcalRecHitsEE"});
0180 desc.add<edm::InputTag>("seedsTag", {"ecalDrivenElectronSeeds"});
0181 desc.add<edm::InputTag>("beamSpotTag", {"offlineBeamSpot"});
0182 desc.add<edm::InputTag>("egmPFCandidatesTag", {"particleFlowEGamma"});
0183
0184
0185 desc.add<bool>("useDefaultEnergyCorrection", true);
0186 desc.add<bool>("useCombinationRegression", false);
0187 desc.add<bool>("ecalDrivenEcalEnergyFromClassBasedParameterization", true);
0188 desc.add<bool>("ecalDrivenEcalErrorFromClassBasedParameterization", true);
0189 desc.add<bool>("applyPreselection", false);
0190 desc.add<bool>("useEcalRegression", false);
0191 desc.add<bool>("applyAmbResolution", false);
0192 desc.add<bool>("ignoreNotPreselected", true);
0193 desc.add<bool>("useGsfPfRecTracks", true);
0194 desc.add<bool>("pureTrackerDrivenEcalErrorFromSimpleParameterization", true);
0195 desc.add<unsigned int>("ambSortingStrategy", 1);
0196 desc.add<unsigned int>("ambClustersOverlapStrategy", 1);
0197 desc.add<bool>("fillConvVtxFitProb", true);
0198 desc.add<bool>("resetMvaValuesUsingPFCandidates", false);
0199
0200
0201 desc.add<std::vector<std::string>>("recHitFlagsToBeExcludedBarrel");
0202 desc.add<std::vector<std::string>>("recHitFlagsToBeExcludedEndcaps");
0203 desc.add<std::vector<std::string>>("recHitSeverityToBeExcludedBarrel");
0204 desc.add<std::vector<std::string>>("recHitSeverityToBeExcludedEndcaps");
0205
0206
0207 desc.add<bool>("checkHcalStatus", true);
0208 desc.add<edm::InputTag>("hbheRecHits", edm::InputTag("hbhereco"));
0209 desc.add<std::vector<double>>("recHitEThresholdHB", {0., 0., 0., 0.});
0210 desc.add<std::vector<double>>("recHitEThresholdHE", {0., 0., 0., 0., 0., 0., 0.});
0211 desc.add<int>("maxHcalRecHitSeverity", 999999);
0212 desc.add<bool>("hcalRun2EffDepth", false);
0213
0214
0215 desc.add("trkIsol03Cfg", EleTkIsolFromCands::pSetDescript());
0216 desc.add("trkIsol04Cfg", EleTkIsolFromCands::pSetDescript());
0217 desc.add("trkIsolHEEP03Cfg", EleTkIsolFromCands::pSetDescript());
0218 desc.add("trkIsolHEEP04Cfg", EleTkIsolFromCands::pSetDescript());
0219 desc.add<bool>("useNumCrystals", true);
0220 desc.add<double>("etMinBarrel", 0.0);
0221 desc.add<double>("etMinEndcaps", 0.11);
0222 desc.add<double>("etMinHcal", 0.0);
0223 desc.add<double>("eMinBarrel", 0.095);
0224 desc.add<double>("eMinEndcaps", 0.0);
0225 desc.add<double>("intRadiusEcalBarrel", 3.0);
0226 desc.add<double>("intRadiusEcalEndcaps", 3.0);
0227 desc.add<double>("intRadiusHcal", 0.15);
0228 desc.add<double>("jurassicWidth", 1.5);
0229 desc.add<bool>("vetoClustered", false);
0230
0231
0232 desc.add<bool>("ctfTracksCheck", true);
0233 desc.add<edm::InputTag>("ctfTracksTag", {"generalTracks"});
0234
0235 desc.add<double>("MaxElePtForOnlyMVA", 50.0);
0236 desc.add<double>("PreSelectMVA", -0.1);
0237
0238 {
0239 edm::ParameterSetDescription psd0;
0240 psd0.add<double>("minSCEtBarrel", 4.0);
0241 psd0.add<double>("minSCEtEndcaps", 4.0);
0242 psd0.add<double>("minEOverPBarrel", 0.0);
0243 psd0.add<double>("minEOverPEndcaps", 0.0);
0244 psd0.add<double>("maxEOverPBarrel", 999999999.0);
0245 psd0.add<double>("maxEOverPEndcaps", 999999999.0);
0246 psd0.add<double>("maxDeltaEtaBarrel", 0.02);
0247 psd0.add<double>("maxDeltaEtaEndcaps", 0.02);
0248 psd0.add<double>("maxDeltaPhiBarrel", 0.15);
0249 psd0.add<double>("maxDeltaPhiEndcaps", 0.15);
0250 psd0.add<double>("hOverEConeSize", 0.15);
0251 psd0.add<double>("maxHOverEBarrelCone", 0.15);
0252 psd0.add<double>("maxHOverEEndcapsCone", 0.15);
0253 psd0.add<double>("maxHBarrelCone", 0.0);
0254 psd0.add<double>("maxHEndcapsCone", 0.0);
0255 psd0.add<double>("maxHOverEBarrelBc", 0.15);
0256 psd0.add<double>("maxHOverEEndcapsBc", 0.15);
0257 psd0.add<double>("maxHBarrelBc", 0.0);
0258 psd0.add<double>("maxHEndcapsBc", 0.0);
0259 psd0.add<double>("maxSigmaIetaIetaBarrel", 999999999.0);
0260 psd0.add<double>("maxSigmaIetaIetaEndcaps", 999999999.0);
0261 psd0.add<double>("maxFbremBarrel", 999999999.0);
0262 psd0.add<double>("maxFbremEndcaps", 999999999.0);
0263 psd0.add<bool>("isBarrel", false);
0264 psd0.add<bool>("isEndcaps", false);
0265 psd0.add<bool>("isFiducial", false);
0266 psd0.add<bool>("seedFromTEC", true);
0267 psd0.add<double>("maxTIP", 999999999.0);
0268 psd0.add<double>("multThresEB", EgammaLocalCovParamDefaults::kMultThresEB);
0269 psd0.add<double>("multThresEE", EgammaLocalCovParamDefaults::kMultThresEE);
0270
0271 desc.add<edm::ParameterSetDescription>("preselection", psd0);
0272 }
0273
0274
0275 desc.add<std::string>("crackCorrectionFunction", "EcalClusterCrackCorrection");
0276
0277 desc.add<bool>("ecalWeightsFromDB", true);
0278 desc.add<std::vector<std::string>>("ecalRefinedRegressionWeightFiles", {})
0279 ->setComment("if not from DB. Otherwise, keep empty");
0280 desc.add<bool>("combinationWeightsFromDB", true);
0281 desc.add<std::vector<std::string>>("combinationRegressionWeightFile", {})
0282 ->setComment("if not from DB. Otherwise, keep empty");
0283
0284
0285 desc.add<std::vector<std::string>>("ecalRefinedRegressionWeightLabels", {});
0286 desc.add<std::vector<std::string>>("combinationRegressionWeightLabels", {});
0287
0288 desc.add<std::vector<std::string>>(
0289 "ElecMVAFilesString",
0290 {
0291 "RecoEgamma/ElectronIdentification/data/TMVA_Category_BDTSimpleCat_10_17Feb2011.weights.xml",
0292 "RecoEgamma/ElectronIdentification/data/TMVA_Category_BDTSimpleCat_12_17Feb2011.weights.xml",
0293 "RecoEgamma/ElectronIdentification/data/TMVA_Category_BDTSimpleCat_20_17Feb2011.weights.xml",
0294 "RecoEgamma/ElectronIdentification/data/TMVA_Category_BDTSimpleCat_22_17Feb2011.weights.xml",
0295 });
0296 desc.add<std::vector<std::string>>(
0297 "SoftElecMVAFilesString",
0298 {
0299 "RecoEgamma/ElectronIdentification/data/TMVA_BDTSoftElectrons_7Feb2014.weights.xml",
0300 });
0301
0302 {
0303 edm::ParameterSetDescription psd1;
0304 psd1.add<bool>("enabled", false);
0305 psd1.add<double>("extetaboundary", 2.65);
0306 psd1.add<std::string>("inputTensorName", "FirstLayer_input");
0307 psd1.add<std::string>("outputTensorName", "sequential/FinalLayer/Softmax");
0308
0309 psd1.add<std::vector<std::string>>(
0310 "modelsFiles",
0311 {"RecoEgamma/ElectronIdentification/data/Ele_PFID_dnn/Run3Summer21_120X/lowpT/lowpT_modelDNN.pb",
0312 "RecoEgamma/ElectronIdentification/data/Ele_PFID_dnn/Run3Summer21_120X/highpTEB/highpTEB_modelDNN.pb",
0313 "RecoEgamma/ElectronIdentification/data/Ele_PFID_dnn/Run3Summer21_120X/highpTEE/highpTEE_modelDNN.pb",
0314 "RecoEgamma/ElectronIdentification/data/Ele_PFID_dnn/Run3Winter22_122X/exteta1/modelDNN.pb",
0315 "RecoEgamma/ElectronIdentification/data/Ele_PFID_dnn/Run3Winter22_122X/exteta2/modelDNN.pb"});
0316 psd1.add<std::vector<std::string>>(
0317 "scalersFiles",
0318 {"RecoEgamma/ElectronIdentification/data/Ele_PFID_dnn/Run3Summer21_120X/lowpT/lowpT_scaler.txt",
0319 "RecoEgamma/ElectronIdentification/data/Ele_PFID_dnn/Run3Summer21_120X/highpTEB/highpTEB_scaler.txt",
0320 "RecoEgamma/ElectronIdentification/data/Ele_PFID_dnn/Run3Summer21_120X/highpTEE/highpTEE_scaler.txt",
0321 "RecoEgamma/ElectronIdentification/data/Ele_PFID_dnn/Run3Winter22_122X/exteta1/scaler.txt",
0322 "RecoEgamma/ElectronIdentification/data/Ele_PFID_dnn/Run3Winter22_122X/exteta2/scaler.txt"});
0323 psd1.add<std::vector<unsigned int>>("outputDim",
0324 {5, 5, 5, 5, 3});
0325
0326 psd1.add<bool>("useEBModelInGap", true);
0327
0328 desc.add<edm::ParameterSetDescription>("EleDNNPFid", psd1);
0329 }
0330
0331
0332
0333 {
0334 edm::ParameterSetDescription psd0;
0335 psd0.add<edm::InputTag>("pfClusterProducer", edm::InputTag("particleFlowClusterECAL"));
0336 psd0.add<double>("drMax", 0.3);
0337 psd0.add<double>("drVetoBarrel", 0.0);
0338 psd0.add<double>("drVetoEndcap", 0.0);
0339 psd0.add<double>("etaStripBarrel", 0.0);
0340 psd0.add<double>("etaStripEndcap", 0.0);
0341 psd0.add<double>("energyBarrel", 0.0);
0342 psd0.add<double>("energyEndcap", 0.0);
0343 desc.add<edm::ParameterSetDescription>("pfECALClusIsolCfg", psd0);
0344 }
0345
0346
0347 {
0348 edm::ParameterSetDescription psd0;
0349 psd0.add<edm::InputTag>("pfClusterProducerHCAL", edm::InputTag("particleFlowClusterHCAL"));
0350 psd0.add<edm::InputTag>("pfClusterProducerHFEM", edm::InputTag(""));
0351 psd0.add<edm::InputTag>("pfClusterProducerHFHAD", edm::InputTag(""));
0352 psd0.add<bool>("useHF", false);
0353 psd0.add<double>("drMax", 0.3);
0354 psd0.add<double>("drVetoBarrel", 0.0);
0355 psd0.add<double>("drVetoEndcap", 0.0);
0356 psd0.add<double>("etaStripBarrel", 0.0);
0357 psd0.add<double>("etaStripEndcap", 0.0);
0358 psd0.add<double>("energyBarrel", 0.0);
0359 psd0.add<double>("energyEndcap", 0.0);
0360 psd0.add<bool>("useEt", true);
0361 desc.add<edm::ParameterSetDescription>("pfHCALClusIsolCfg", psd0);
0362 }
0363
0364 descriptions.add("gsfElectronProducerDefault", desc);
0365 }
0366
0367 namespace {
0368 GsfElectronAlgo::CutsConfiguration makeCutsConfiguration(edm::ParameterSet const& pset) {
0369 return GsfElectronAlgo::CutsConfiguration{
0370 .minSCEtBarrel = pset.getParameter<double>("minSCEtBarrel"),
0371 .minSCEtEndcaps = pset.getParameter<double>("minSCEtEndcaps"),
0372 .maxEOverPBarrel = pset.getParameter<double>("maxEOverPBarrel"),
0373 .maxEOverPEndcaps = pset.getParameter<double>("maxEOverPEndcaps"),
0374 .minEOverPBarrel = pset.getParameter<double>("minEOverPBarrel"),
0375 .minEOverPEndcaps = pset.getParameter<double>("minEOverPEndcaps"),
0376 .maxHOverEBarrelCone = pset.getParameter<double>("maxHOverEBarrelCone"),
0377 .maxHOverEEndcapsCone = pset.getParameter<double>("maxHOverEEndcapsCone"),
0378 .maxHBarrelCone = pset.getParameter<double>("maxHBarrelCone"),
0379 .maxHEndcapsCone = pset.getParameter<double>("maxHEndcapsCone"),
0380 .maxHOverEBarrelBc = pset.getParameter<double>("maxHOverEBarrelBc"),
0381 .maxHOverEEndcapsBc = pset.getParameter<double>("maxHOverEEndcapsBc"),
0382 .maxHBarrelBc = pset.getParameter<double>("maxHBarrelBc"),
0383 .maxHEndcapsBc = pset.getParameter<double>("maxHEndcapsBc"),
0384 .maxDeltaEtaBarrel = pset.getParameter<double>("maxDeltaEtaBarrel"),
0385 .maxDeltaEtaEndcaps = pset.getParameter<double>("maxDeltaEtaEndcaps"),
0386 .maxDeltaPhiBarrel = pset.getParameter<double>("maxDeltaPhiBarrel"),
0387 .maxDeltaPhiEndcaps = pset.getParameter<double>("maxDeltaPhiEndcaps"),
0388 .maxSigmaIetaIetaBarrel = pset.getParameter<double>("maxSigmaIetaIetaBarrel"),
0389 .maxSigmaIetaIetaEndcaps = pset.getParameter<double>("maxSigmaIetaIetaEndcaps"),
0390 .maxFbremBarrel = pset.getParameter<double>("maxFbremBarrel"),
0391 .maxFbremEndcaps = pset.getParameter<double>("maxFbremEndcaps"),
0392 .isBarrel = pset.getParameter<bool>("isBarrel"),
0393 .isEndcaps = pset.getParameter<bool>("isEndcaps"),
0394 .isFiducial = pset.getParameter<bool>("isFiducial"),
0395 .maxTIP = pset.getParameter<double>("maxTIP"),
0396 .seedFromTEC = pset.getParameter<bool>("seedFromTEC"),
0397 .multThresEB = pset.getParameter<double>("multThresEB"),
0398 .multThresEE = pset.getParameter<double>("multThresEE"),
0399 };
0400 }
0401 };
0402
0403 GsfElectronProducer::GsfElectronProducer(const edm::ParameterSet& cfg, const GsfElectronAlgo::HeavyObjectCache* gcache)
0404 : cutsCfg_{makeCutsConfiguration(cfg.getParameter<edm::ParameterSet>("preselection"))},
0405 ecalSeedingParametersChecked_(false),
0406 electronPutToken_(produces<GsfElectronCollection>()),
0407 gsfPfRecTracksTag_(consumes(cfg.getParameter<edm::InputTag>("gsfPfRecTracksTag"))),
0408 useGsfPfRecTracks_(cfg.getParameter<bool>("useGsfPfRecTracks")),
0409 resetMvaValuesUsingPFCandidates_(cfg.getParameter<bool>("resetMvaValuesUsingPFCandidates")) {
0410 if (resetMvaValuesUsingPFCandidates_) {
0411 egmPFCandidateCollection_ = consumes(cfg.getParameter<edm::InputTag>("egmPFCandidatesTag"));
0412 }
0413
0414 inputCfg_.gsfElectronCores = consumes(cfg.getParameter<edm::InputTag>("gsfElectronCoresTag"));
0415 inputCfg_.hbheRecHitsTag = consumes(cfg.getParameter<edm::InputTag>("hbheRecHits"));
0416 inputCfg_.barrelRecHitCollection = consumes(cfg.getParameter<edm::InputTag>("barrelRecHitCollectionTag"));
0417 inputCfg_.endcapRecHitCollection = consumes(cfg.getParameter<edm::InputTag>("endcapRecHitCollectionTag"));
0418 inputCfg_.ctfTracks = consumes(cfg.getParameter<edm::InputTag>("ctfTracksTag"));
0419
0420 inputCfg_.seedsTag = consumes(cfg.getParameter<edm::InputTag>("seedsTag"));
0421 inputCfg_.beamSpotTag = consumes(cfg.getParameter<edm::InputTag>("beamSpotTag"));
0422 inputCfg_.vtxCollectionTag = consumes(cfg.getParameter<edm::InputTag>("vtxTag"));
0423 if (cfg.getParameter<bool>("fillConvVtxFitProb"))
0424 inputCfg_.conversions = consumes(cfg.getParameter<edm::InputTag>("conversionsTag"));
0425
0426
0427 const edm::ParameterSet& pfECALClusIsolCfg = cfg.getParameter<edm::ParameterSet>("pfECALClusIsolCfg");
0428 const edm::ParameterSet& pfHCALClusIsolCfg = cfg.getParameter<edm::ParameterSet>("pfHCALClusIsolCfg");
0429 inputCfg_.pfClusterProducer =
0430 consumes<reco::PFClusterCollection>(pfECALClusIsolCfg.getParameter<edm::InputTag>("pfClusterProducer"));
0431 inputCfg_.pfClusterProducerHCAL = consumes(pfHCALClusIsolCfg.getParameter<edm::InputTag>("pfClusterProducerHCAL"));
0432 inputCfg_.pfClusterProducerHFEM = consumes(pfHCALClusIsolCfg.getParameter<edm::InputTag>("pfClusterProducerHFEM"));
0433 inputCfg_.pfClusterProducerHFHAD = consumes(pfHCALClusIsolCfg.getParameter<edm::InputTag>("pfClusterProducerHFHAD"));
0434
0435
0436 const auto& pset_dnn = cfg.getParameter<edm::ParameterSet>("EleDNNPFid");
0437 dnnPFidEnabled_ = pset_dnn.getParameter<bool>("enabled");
0438 extetaboundary_ = pset_dnn.getParameter<double>("extetaboundary");
0439
0440 strategyCfg_.useDefaultEnergyCorrection = cfg.getParameter<bool>("useDefaultEnergyCorrection");
0441
0442 strategyCfg_.applyPreselection = cfg.getParameter<bool>("applyPreselection");
0443 strategyCfg_.ecalDrivenEcalEnergyFromClassBasedParameterization =
0444 cfg.getParameter<bool>("ecalDrivenEcalEnergyFromClassBasedParameterization");
0445 strategyCfg_.ecalDrivenEcalErrorFromClassBasedParameterization =
0446 cfg.getParameter<bool>("ecalDrivenEcalErrorFromClassBasedParameterization");
0447 strategyCfg_.pureTrackerDrivenEcalErrorFromSimpleParameterization =
0448 cfg.getParameter<bool>("pureTrackerDrivenEcalErrorFromSimpleParameterization");
0449 strategyCfg_.applyAmbResolution = cfg.getParameter<bool>("applyAmbResolution");
0450 strategyCfg_.ignoreNotPreselected = cfg.getParameter<bool>("ignoreNotPreselected");
0451 strategyCfg_.ambSortingStrategy = cfg.getParameter<unsigned>("ambSortingStrategy");
0452 strategyCfg_.ambClustersOverlapStrategy = cfg.getParameter<unsigned>("ambClustersOverlapStrategy");
0453 strategyCfg_.ctfTracksCheck = cfg.getParameter<bool>("ctfTracksCheck");
0454 strategyCfg_.PreSelectMVA = cfg.getParameter<double>("PreSelectMVA");
0455 strategyCfg_.MaxElePtForOnlyMVA = cfg.getParameter<double>("MaxElePtForOnlyMVA");
0456 strategyCfg_.useEcalRegression = cfg.getParameter<bool>("useEcalRegression");
0457 strategyCfg_.useCombinationRegression = cfg.getParameter<bool>("useCombinationRegression");
0458 strategyCfg_.fillConvVtxFitProb = cfg.getParameter<bool>("fillConvVtxFitProb");
0459 strategyCfg_.computePfClusterIso = dnnPFidEnabled_;
0460
0461
0462 auto const& psetPreselection = cfg.getParameter<edm::ParameterSet>("preselection");
0463 hcalCfg_.hOverEConeSize = psetPreselection.getParameter<double>("hOverEConeSize");
0464 if (hcalCfg_.hOverEConeSize > 0) {
0465 hcalCfg_.onlyBehindCluster = false;
0466 hcalCfg_.checkHcalStatus = cfg.getParameter<bool>("checkHcalStatus");
0467
0468
0469 hcalCfg_.hbheRecHits = consumes<HBHERecHitCollection>(cfg.getParameter<edm::InputTag>("hbheRecHits"));
0470
0471 hcalCfg_.eThresHB = cfg.getParameter<EgammaHcalIsolation::arrayHB>("recHitEThresholdHB");
0472 hcalCfg_.maxSeverityHB = cfg.getParameter<int>("maxHcalRecHitSeverity");
0473 hcalCfg_.eThresHE = cfg.getParameter<EgammaHcalIsolation::arrayHE>("recHitEThresholdHE");
0474 hcalCfg_.maxSeverityHE = hcalCfg_.maxSeverityHB;
0475 }
0476
0477 hcalCfgBc_.hOverEConeSize = 0.;
0478 hcalCfgBc_.onlyBehindCluster = true;
0479 hcalCfgBc_.checkHcalStatus = cfg.getParameter<bool>("checkHcalStatus");
0480
0481
0482 hcalCfgBc_.hbheRecHits = consumes<HBHERecHitCollection>(cfg.getParameter<edm::InputTag>("hbheRecHits"));
0483
0484 hcalCfgBc_.eThresHB = cfg.getParameter<EgammaHcalIsolation::arrayHB>("recHitEThresholdHB");
0485 hcalCfgBc_.maxSeverityHB = cfg.getParameter<int>("maxHcalRecHitSeverity");
0486 hcalCfgBc_.eThresHE = cfg.getParameter<EgammaHcalIsolation::arrayHE>("recHitEThresholdHE");
0487 hcalCfgBc_.maxSeverityHE = hcalCfgBc_.maxSeverityHB;
0488
0489 hcalRun2EffDepth_ = cfg.getParameter<bool>("hcalRun2EffDepth");
0490
0491
0492 GsfElectronAlgo::EcalRecHitsConfiguration recHitsCfg;
0493 auto const& flagnamesbarrel = cfg.getParameter<std::vector<std::string>>("recHitFlagsToBeExcludedBarrel");
0494 recHitsCfg.recHitFlagsToBeExcludedBarrel = StringToEnumValue<EcalRecHit::Flags>(flagnamesbarrel);
0495 auto const& flagnamesendcaps = cfg.getParameter<std::vector<std::string>>("recHitFlagsToBeExcludedEndcaps");
0496 recHitsCfg.recHitFlagsToBeExcludedEndcaps = StringToEnumValue<EcalRecHit::Flags>(flagnamesendcaps);
0497 auto const& severitynamesbarrel = cfg.getParameter<std::vector<std::string>>("recHitSeverityToBeExcludedBarrel");
0498 recHitsCfg.recHitSeverityToBeExcludedBarrel =
0499 StringToEnumValue<EcalSeverityLevel::SeverityLevel>(severitynamesbarrel);
0500 auto const& severitynamesendcaps = cfg.getParameter<std::vector<std::string>>("recHitSeverityToBeExcludedEndcaps");
0501 recHitsCfg.recHitSeverityToBeExcludedEndcaps =
0502 StringToEnumValue<EcalSeverityLevel::SeverityLevel>(severitynamesendcaps);
0503
0504
0505 const GsfElectronAlgo::IsolationConfiguration isoCfg{
0506 .intRadiusHcal = cfg.getParameter<double>("intRadiusHcal"),
0507 .etMinHcal = cfg.getParameter<double>("etMinHcal"),
0508 .intRadiusEcalBarrel = cfg.getParameter<double>("intRadiusEcalBarrel"),
0509 .intRadiusEcalEndcaps = cfg.getParameter<double>("intRadiusEcalEndcaps"),
0510 .jurassicWidth = cfg.getParameter<double>("jurassicWidth"),
0511 .etMinBarrel = cfg.getParameter<double>("etMinBarrel"),
0512 .eMinBarrel = cfg.getParameter<double>("eMinBarrel"),
0513 .etMinEndcaps = cfg.getParameter<double>("etMinEndcaps"),
0514 .eMinEndcaps = cfg.getParameter<double>("eMinEndcaps"),
0515 .vetoClustered = cfg.getParameter<bool>("vetoClustered"),
0516 .useNumCrystals = cfg.getParameter<bool>("useNumCrystals")};
0517
0518
0519 const GsfElectronAlgo::PFClusterIsolationConfiguration pfisoCfg{
0520 .ecaldrMax = pfECALClusIsolCfg.getParameter<double>("drMax"),
0521 .ecaldrVetoBarrel = pfECALClusIsolCfg.getParameter<double>("drVetoBarrel"),
0522 .ecaldrVetoEndcap = pfECALClusIsolCfg.getParameter<double>("drVetoEndcap"),
0523 .ecaletaStripBarrel = pfECALClusIsolCfg.getParameter<double>("etaStripBarrel"),
0524 .ecaletaStripEndcap = pfECALClusIsolCfg.getParameter<double>("etaStripEndcap"),
0525 .ecalenergyBarrel = pfECALClusIsolCfg.getParameter<double>("energyBarrel"),
0526 .ecalenergyEndcap = pfECALClusIsolCfg.getParameter<double>("energyEndcap"),
0527 .useHF = pfHCALClusIsolCfg.getParameter<bool>("useHF"),
0528 .hcaldrMax = pfHCALClusIsolCfg.getParameter<double>("drMax"),
0529 .hcaldrVetoBarrel = pfHCALClusIsolCfg.getParameter<double>("drVetoBarrel"),
0530 .hcaldrVetoEndcap = pfHCALClusIsolCfg.getParameter<double>("drVetoEndcap"),
0531 .hcaletaStripBarrel = pfHCALClusIsolCfg.getParameter<double>("etaStripBarrel"),
0532 .hcaletaStripEndcap = pfHCALClusIsolCfg.getParameter<double>("etaStripEndcap"),
0533 .hcalenergyBarrel = pfHCALClusIsolCfg.getParameter<double>("energyBarrel"),
0534 .hcalenergyEndcap = pfHCALClusIsolCfg.getParameter<double>("energyEndcap"),
0535 .hcaluseEt = pfHCALClusIsolCfg.getParameter<bool>("useEt")};
0536
0537 const RegressionHelper::Configuration regressionCfg{
0538 .ecalRegressionWeightLabels = cfg.getParameter<std::vector<std::string>>("ecalRefinedRegressionWeightLabels"),
0539 .ecalWeightsFromDB = cfg.getParameter<bool>("ecalWeightsFromDB"),
0540 .ecalRegressionWeightFiles = cfg.getParameter<std::vector<std::string>>("ecalRefinedRegressionWeightFiles"),
0541 .combinationRegressionWeightLabels =
0542 cfg.getParameter<std::vector<std::string>>("combinationRegressionWeightLabels"),
0543 .combinationWeightsFromDB = cfg.getParameter<bool>("combinationWeightsFromDB"),
0544 .combinationRegressionWeightFiles =
0545 cfg.getParameter<std::vector<std::string>>("combinationRegressionWeightFile")};
0546
0547
0548 algo_ = std::make_unique<GsfElectronAlgo>(
0549 inputCfg_,
0550 strategyCfg_,
0551 cutsCfg_,
0552 hcalCfg_,
0553 hcalCfgBc_,
0554 isoCfg,
0555 pfisoCfg,
0556 recHitsCfg,
0557 EcalClusterFunctionFactory::get()->create(
0558 cfg.getParameter<std::string>("crackCorrectionFunction"), cfg, consumesCollector()),
0559 regressionCfg,
0560 cfg.getParameter<edm::ParameterSet>("trkIsol03Cfg"),
0561 cfg.getParameter<edm::ParameterSet>("trkIsol04Cfg"),
0562 cfg.getParameter<edm::ParameterSet>("trkIsolHEEP03Cfg"),
0563 cfg.getParameter<edm::ParameterSet>("trkIsolHEEP04Cfg"),
0564 consumesCollector());
0565
0566 if (dnnPFidEnabled_) {
0567 tfSessions_ = gcache->iElectronDNNEstimator->getSessions();
0568 }
0569 }
0570
0571 void GsfElectronProducer::endStream() {
0572 for (auto session : tfSessions_) {
0573 tensorflow::closeSession(session);
0574 }
0575 }
0576
0577 void GsfElectronProducer::checkEcalSeedingParameters(edm::ParameterSet const& pset) {
0578 if (!pset.exists("SeedConfiguration")) {
0579 return;
0580 }
0581 edm::ParameterSet seedConfiguration = pset.getParameter<edm::ParameterSet>("SeedConfiguration");
0582
0583 if (seedConfiguration.getParameter<bool>("applyHOverECut")) {
0584 if ((hcalCfg_.hOverEConeSize != 0) &&
0585 (hcalCfg_.hOverEConeSize != seedConfiguration.getParameter<double>("hOverEConeSize"))) {
0586 edm::LogWarning("GsfElectronAlgo|InconsistentParameters")
0587 << "The H/E cone size (" << hcalCfg_.hOverEConeSize << ") is different from ecal seeding ("
0588 << seedConfiguration.getParameter<double>("hOverEConeSize") << ").";
0589 }
0590 if (cutsCfg_.maxHOverEBarrelCone < seedConfiguration.getParameter<double>("maxHOverEBarrel")) {
0591 edm::LogWarning("GsfElectronAlgo|InconsistentParameters")
0592 << "The max barrel cone H/E is lower than during ecal seeding.";
0593 }
0594 if (cutsCfg_.maxHOverEEndcapsCone < seedConfiguration.getParameter<double>("maxHOverEEndcaps")) {
0595 edm::LogWarning("GsfElectronAlgo|InconsistentParameters")
0596 << "The max endcaps cone H/E is lower than during ecal seeding.";
0597 }
0598 }
0599
0600 if (cutsCfg_.minSCEtBarrel < seedConfiguration.getParameter<double>("SCEtCut")) {
0601 edm::LogWarning("GsfElectronAlgo|InconsistentParameters")
0602 << "The minimum super-cluster Et in barrel is lower than during ecal seeding.";
0603 }
0604 if (cutsCfg_.minSCEtEndcaps < seedConfiguration.getParameter<double>("SCEtCut")) {
0605 edm::LogWarning("GsfElectronAlgo|InconsistentParameters")
0606 << "The minimum super-cluster Et in endcaps is lower than during ecal seeding.";
0607 }
0608 }
0609
0610
0611
0612
0613
0614 void GsfElectronProducer::setAmbiguityData(reco::GsfElectronCollection& electrons,
0615 edm::Event const& event,
0616 bool ignoreNotPreselected) const {
0617
0618 auto const& beamspot = event.get(inputCfg_.beamSpotTag);
0619 auto gsfPfRecTracks =
0620 useGsfPfRecTracks_ ? event.getHandle(gsfPfRecTracksTag_) : edm::Handle<reco::GsfPFRecTrackCollection>{};
0621 auto const& barrelRecHits = event.get(inputCfg_.barrelRecHitCollection);
0622 auto const& endcapRecHits = event.get(inputCfg_.endcapRecHitCollection);
0623
0624 if (strategyCfg_.ambSortingStrategy == 0) {
0625 std::sort(electrons.begin(), electrons.end(), egamma::isBetterElectron);
0626 } else if (strategyCfg_.ambSortingStrategy == 1) {
0627 std::sort(electrons.begin(), electrons.end(), egamma::isInnermostElectron);
0628 } else {
0629 throw cms::Exception("GsfElectronAlgo|UnknownAmbiguitySortingStrategy")
0630 << "value of strategyCfg_.ambSortingStrategy is : " << strategyCfg_.ambSortingStrategy;
0631 }
0632
0633
0634 for (auto& electron : electrons) {
0635 electron.clearAmbiguousGsfTracks();
0636 electron.setAmbiguous(false);
0637 }
0638
0639
0640 if (useGsfPfRecTracks_) {
0641 for (auto& e1 : electrons) {
0642 bool found = false;
0643 for (auto const& gsfPfRecTrack : *gsfPfRecTracks) {
0644 if (gsfPfRecTrack.gsfTrackRef() == e1.gsfTrack()) {
0645 if (found) {
0646 edm::LogWarning("GsfElectronAlgo") << "associated gsfPfRecTrack already found";
0647 } else {
0648 found = true;
0649 for (auto const& duplicate : gsfPfRecTrack.convBremGsfPFRecTrackRef()) {
0650 e1.addAmbiguousGsfTrack(duplicate->gsfTrackRef());
0651 }
0652 }
0653 }
0654 }
0655 }
0656 }
0657
0658 else {
0659 for (auto e1 = electrons.begin(); e1 != electrons.end(); ++e1) {
0660 if (e1->ambiguous())
0661 continue;
0662 if (ignoreNotPreselected && !isPreselected(*e1))
0663 continue;
0664
0665 SuperClusterRef scRef1 = e1->superCluster();
0666 CaloClusterPtr eleClu1 = e1->electronCluster();
0667 LogDebug("GsfElectronAlgo") << "Blessing electron with E/P " << e1->eSuperClusterOverP() << ", cluster "
0668 << scRef1.get() << " & track " << e1->gsfTrack().get();
0669
0670 for (auto e2 = e1 + 1; e2 != electrons.end(); ++e2) {
0671 if (e2->ambiguous())
0672 continue;
0673 if (ignoreNotPreselected && !isPreselected(*e2))
0674 continue;
0675
0676 SuperClusterRef scRef2 = e2->superCluster();
0677 CaloClusterPtr eleClu2 = e2->electronCluster();
0678
0679
0680 bool sameCluster = false;
0681 if (strategyCfg_.ambClustersOverlapStrategy == 0) {
0682 sameCluster = (scRef1 == scRef2);
0683 } else if (strategyCfg_.ambClustersOverlapStrategy == 1) {
0684 float eMin = 1.;
0685 float threshold = eMin * cosh(EleRelPoint(scRef1->position(), beamspot.position()).eta());
0686 using egamma::sharedEnergy;
0687 sameCluster = ((sharedEnergy(*eleClu1, *eleClu2, barrelRecHits, endcapRecHits) >= threshold) ||
0688 (sharedEnergy(*scRef1->seed(), *eleClu2, barrelRecHits, endcapRecHits) >= threshold) ||
0689 (sharedEnergy(*eleClu1, *scRef2->seed(), barrelRecHits, endcapRecHits) >= threshold) ||
0690 (sharedEnergy(*scRef1->seed(), *scRef2->seed(), barrelRecHits, endcapRecHits) >= threshold));
0691 } else {
0692 throw cms::Exception("GsfElectronAlgo|UnknownAmbiguityClustersOverlapStrategy")
0693 << "value of strategyCfg_.ambClustersOverlapStrategy is : " << strategyCfg_.ambClustersOverlapStrategy;
0694 }
0695
0696
0697 if (sameCluster) {
0698 LogDebug("GsfElectronAlgo") << "Discarding electron with E/P " << e2->eSuperClusterOverP() << ", cluster "
0699 << scRef2.get() << " and track " << e2->gsfTrack().get();
0700 e1->addAmbiguousGsfTrack(e2->gsfTrack());
0701 e2->setAmbiguous(true);
0702 } else if (e1->gsfTrack() == e2->gsfTrack()) {
0703 edm::LogWarning("GsfElectronAlgo") << "Forgetting electron with E/P " << e2->eSuperClusterOverP()
0704 << ", cluster " << scRef2.get() << " and track " << e2->gsfTrack().get();
0705 e2->setAmbiguous(true);
0706 }
0707 }
0708 }
0709 }
0710 }
0711
0712 bool GsfElectronProducer::isPreselected(GsfElectron const& ele) const {
0713 bool passCutBased = ele.passingCutBasedPreselection();
0714 bool passPF = ele.passingPflowPreselection();
0715
0716
0717 bool passmva = ele.passingMvaPreselection();
0718 if (!ele.ecalDrivenSeed()) {
0719 if (ele.pt() > strategyCfg_.MaxElePtForOnlyMVA)
0720 return passmva && passCutBased;
0721 else
0722 return passmva;
0723 } else {
0724 return passCutBased || passPF || passmva;
0725 }
0726 }
0727
0728
0729 void GsfElectronProducer::produce(edm::Event& event, const edm::EventSetup& setup) {
0730
0731 if (!ecalSeedingParametersChecked_) {
0732 ecalSeedingParametersChecked_ = true;
0733 auto seeds = event.getHandle(inputCfg_.seedsTag);
0734 if (!seeds.isValid()) {
0735 edm::LogWarning("GsfElectronAlgo|UnreachableSeedsProvenance")
0736 << "Cannot check consistency of parameters with ecal seeding ones,"
0737 << " because the original collection of seeds is not any more available.";
0738 } else {
0739 checkEcalSeedingParameters(edm::parameterSet(seeds.provenance()->stable(), event.processHistory()));
0740 }
0741 }
0742
0743 auto electrons = algo_->completeElectrons(event, setup, globalCache());
0744 if (resetMvaValuesUsingPFCandidates_) {
0745 const auto gsfMVAInputMap = matchWithPFCandidates(event.get(egmPFCandidateCollection_));
0746 for (auto& el : electrons) {
0747 el.setMvaInput(gsfMVAInputMap.find(el.gsfTrack())->second);
0748 }
0749 setMVAOutputs(
0750 electrons, globalCache(), event.get(inputCfg_.vtxCollectionTag), dnnPFidEnabled_, extetaboundary_, tfSessions_);
0751 }
0752
0753
0754 logElectrons(electrons, event, "GsfElectronAlgo Info (before preselection)");
0755
0756 if (strategyCfg_.applyPreselection) {
0757 electrons.erase(
0758 std::remove_if(electrons.begin(), electrons.end(), [this](auto const& ele) { return !isPreselected(ele); }),
0759 electrons.end());
0760 logElectrons(electrons, event, "GsfElectronAlgo Info (after preselection)");
0761 }
0762
0763 setAmbiguityData(electrons, event, strategyCfg_.ignoreNotPreselected);
0764 if (strategyCfg_.applyAmbResolution) {
0765 electrons.erase(std::remove_if(electrons.begin(), electrons.end(), std::mem_fn(&reco::GsfElectron::ambiguous)),
0766 electrons.end());
0767 logElectrons(electrons, event, "GsfElectronAlgo Info (after amb. solving)");
0768 }
0769
0770 if (hcalRun2EffDepth_) {
0771 for (auto& ele : electrons)
0772 ele.hcalToRun2EffDepth();
0773 }
0774
0775 event.emplace(electronPutToken_, std::move(electrons));
0776 }
0777
0778 #include "FWCore/Framework/interface/MakerMacros.h"
0779 DEFINE_FWK_MODULE(GsfElectronProducer);