File indexing completed on 2022-06-16 01:50:15
0001
0002
0003
0004
0005
0006
0007
0008 #include "CommonTools/Utils/interface/StringToEnumValue.h"
0009 #include "CondFormats/DataRecord/interface/HcalChannelQualityRcd.h"
0010 #include "CondFormats/EcalObjects/interface/EcalFunctionParameters.h"
0011 #include "DataFormats/Common/interface/Handle.h"
0012 #include "DataFormats/Common/interface/ValueMap.h"
0013 #include "DataFormats/EgammaCandidates/interface/Conversion.h"
0014 #include "DataFormats/EgammaCandidates/interface/Photon.h"
0015 #include "DataFormats/EgammaCandidates/interface/PhotonCore.h"
0016 #include "DataFormats/EgammaCandidates/interface/PhotonFwd.h"
0017 #include "DataFormats/EgammaReco/interface/BasicCluster.h"
0018 #include "DataFormats/EgammaReco/interface/BasicClusterShapeAssociation.h"
0019 #include "DataFormats/EgammaReco/interface/ClusterShape.h"
0020 #include "DataFormats/EgammaReco/interface/SuperCluster.h"
0021 #include "DataFormats/ParticleFlowCandidate/interface/PFCandidate.h"
0022 #include "DataFormats/ParticleFlowCandidate/interface/PFCandidateEGammaExtra.h"
0023 #include "DataFormats/VertexReco/interface/Vertex.h"
0024 #include "FWCore/Framework/interface/ESHandle.h"
0025 #include "FWCore/Framework/interface/Event.h"
0026 #include "FWCore/Framework/interface/EventSetup.h"
0027 #include "FWCore/Framework/interface/stream/EDProducer.h"
0028 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0029 #include "FWCore/ParameterSet/interface/ParameterSet.h"
0030 #include "FWCore/Utilities/interface/Exception.h"
0031 #include "FWCore/Utilities/interface/isFinite.h"
0032 #include "Geometry/Records/interface/CaloTopologyRecord.h"
0033 #include "Geometry/CaloGeometry/interface/CaloGeometry.h"
0034 #include "Geometry/CaloGeometry/interface/CaloSubdetectorGeometry.h"
0035 #include "Geometry/CaloTopology/interface/CaloTopology.h"
0036 #include "Geometry/Records/interface/CaloGeometryRecord.h"
0037 #include "RecoEcal/EgammaCoreTools/interface/EcalClusterFunctionBaseClass.h"
0038 #include "RecoEcal/EgammaCoreTools/interface/EcalClusterTools.h"
0039 #include "RecoEcal/EgammaCoreTools/interface/EcalTools.h"
0040 #include "RecoEgamma/EgammaPhotonAlgos/interface/PhotonEnergyCorrector.h"
0041 #include "RecoEgamma/PhotonIdentification/interface/PhotonIsolationCalculator.h"
0042 #include "RecoEgamma/PhotonIdentification/interface/PhotonMIPHaloTagger.h"
0043 #include "RecoEgamma/PhotonIdentification/interface/PhotonMVABasedHaloTagger.h"
0044 #include "RecoLocalCalo/EcalRecAlgos/interface/EcalSeverityLevelAlgo.h"
0045 #include "RecoLocalCalo/EcalRecAlgos/interface/EcalSeverityLevelAlgoRcd.h"
0046 #include "CondFormats/EcalObjects/interface/EcalPFRecHitThresholds.h"
0047 #include "CondFormats/DataRecord/interface/EcalPFRecHitThresholdsRcd.h"
0048 #include "RecoEcal/EgammaCoreTools/interface/EgammaLocalCovParamDefaults.h"
0049 #include "RecoEgamma/EgammaElectronAlgos/interface/ElectronHcalHelper.h"
0050 #include "RecoEgamma/PhotonIdentification/interface/PhotonDNNEstimator.h"
0051 #include "PhysicsTools/TensorFlow/interface/TensorFlow.h"
0052 #include "RecoEgamma/EgammaIsolationAlgos/interface/EcalPFClusterIsolation.h"
0053 #include "RecoEgamma/EgammaIsolationAlgos/interface/HcalPFClusterIsolation.h"
0054 #include "CondFormats/GBRForest/interface/GBRForest.h"
0055 #include "CommonTools/MVAUtils/interface/GBRForestTools.h"
0056
0057 class CacheData {
0058 public:
0059 CacheData(const edm::ParameterSet& conf) {
0060
0061 egammaTools::DNNConfiguration config;
0062 const auto& pset_dnn = conf.getParameter<edm::ParameterSet>("PhotonDNNPFid");
0063 const auto dnnEnabled = pset_dnn.getParameter<bool>("enabled");
0064 if (dnnEnabled) {
0065 config.inputTensorName = pset_dnn.getParameter<std::string>("inputTensorName");
0066 config.outputTensorName = pset_dnn.getParameter<std::string>("outputTensorName");
0067 config.modelsFiles = pset_dnn.getParameter<std::vector<std::string>>("modelsFiles");
0068 config.scalersFiles = pset_dnn.getParameter<std::vector<std::string>>("scalersFiles");
0069 config.outputDim = pset_dnn.getParameter<std::vector<unsigned int>>("outputDim");
0070 const auto useEBModelInGap = pset_dnn.getParameter<bool>("useEBModelInGap");
0071 photonDNNEstimator = std::make_unique<PhotonDNNEstimator>(config, useEBModelInGap);
0072 }
0073
0074 const auto runMVABasedHaloTagger = conf.getParameter<bool>("runMVABasedHaloTagger");
0075 edm::ParameterSet mvaBasedHaloVariableSet = conf.getParameter<edm::ParameterSet>("mvaBasedHaloVariableSet");
0076 auto trainingFileName_ = mvaBasedHaloVariableSet.getParameter<edm::FileInPath>(("trainingFileName")).fullPath();
0077 if (runMVABasedHaloTagger) {
0078 haloTaggerGBR = createGBRForest(trainingFileName_);
0079 }
0080 }
0081 std::unique_ptr<const PhotonDNNEstimator> photonDNNEstimator;
0082 std::unique_ptr<const GBRForest> haloTaggerGBR;
0083 };
0084
0085 class GEDPhotonProducer : public edm::stream::EDProducer<edm::GlobalCache<CacheData>> {
0086 public:
0087 GEDPhotonProducer(const edm::ParameterSet& ps, const CacheData* gcache);
0088
0089 void produce(edm::Event& evt, const edm::EventSetup& es) override;
0090
0091 static std::unique_ptr<CacheData> initializeGlobalCache(const edm::ParameterSet&);
0092 static void globalEndJob(const CacheData*){};
0093
0094 void endStream() override;
0095
0096 private:
0097 class RecoStepInfo {
0098 public:
0099 enum FlagBits { kOOT = 0x1, kFinal = 0x2 };
0100 explicit RecoStepInfo(const std::string& recoStep);
0101
0102 bool isOOT() const { return flags_ & kOOT; }
0103 bool isFinal() const { return flags_ & kFinal; }
0104
0105 private:
0106 unsigned int flags_;
0107 };
0108
0109 void fillPhotonCollection(edm::Event& evt,
0110 edm::EventSetup const& es,
0111 const edm::Handle<reco::PhotonCoreCollection>& photonCoreHandle,
0112 const CaloTopology* topology,
0113 const EcalRecHitCollection* ecalBarrelHits,
0114 const EcalRecHitCollection* ecalEndcapHits,
0115 const EcalRecHitCollection* preshowerHits,
0116 const ElectronHcalHelper* hcalHelperCone,
0117 const ElectronHcalHelper* hcalHelperBc,
0118 const reco::VertexCollection& pvVertices,
0119 reco::PhotonCollection& outputCollection,
0120 int& iSC,
0121 EcalPFRecHitThresholds const& thresholds);
0122
0123 void fillPhotonCollection(edm::Event& evt,
0124 edm::EventSetup const& es,
0125 const edm::Handle<reco::PhotonCollection>& photonHandle,
0126 const edm::Handle<reco::PFCandidateCollection> pfCandidateHandle,
0127 const edm::Handle<reco::PFCandidateCollection> pfEGCandidateHandle,
0128 reco::VertexCollection const& pvVertices,
0129 reco::PhotonCollection& outputCollection,
0130 int& iSC,
0131 const edm::Handle<edm::ValueMap<float>>& chargedHadrons,
0132 const edm::Handle<edm::ValueMap<float>>& neutralHadrons,
0133 const edm::Handle<edm::ValueMap<float>>& photons,
0134 const edm::Handle<edm::ValueMap<float>>& chargedHadronsWorstVtx,
0135 const edm::Handle<edm::ValueMap<float>>& chargedHadronsWorstVtxGeomVeto,
0136 const edm::Handle<edm::ValueMap<float>>& chargedHadronsPFPV,
0137 const edm::Handle<edm::ValueMap<float>>& pfEcalClusters,
0138 const edm::Handle<edm::ValueMap<float>>& pfHcalClusters);
0139
0140
0141 std::string photonCollection_;
0142 const edm::InputTag photonProducer_;
0143
0144 edm::EDGetTokenT<reco::PhotonCoreCollection> photonCoreProducerT_;
0145 edm::EDGetTokenT<reco::PhotonCollection> photonProducerT_;
0146 edm::EDGetTokenT<EcalRecHitCollection> barrelEcalHits_;
0147 edm::EDGetTokenT<EcalRecHitCollection> endcapEcalHits_;
0148 edm::EDGetTokenT<EcalRecHitCollection> preshowerHits_;
0149 edm::EDGetTokenT<reco::PFCandidateCollection> pfEgammaCandidates_;
0150 edm::EDGetTokenT<reco::PFCandidateCollection> pfCandidates_;
0151 edm::EDGetTokenT<HBHERecHitCollection> hbheRecHits_;
0152 edm::EDGetTokenT<reco::VertexCollection> vertexProducer_;
0153
0154 edm::EDGetTokenT<edm::ValueMap<std::vector<reco::PFCandidateRef>>> particleBasedIsolationToken;
0155
0156 edm::EDGetTokenT<edm::ValueMap<float>> phoChargedIsolationToken_;
0157 edm::EDGetTokenT<edm::ValueMap<float>> phoNeutralHadronIsolationToken_;
0158 edm::EDGetTokenT<edm::ValueMap<float>> phoPhotonIsolationToken_;
0159 edm::EDGetTokenT<edm::ValueMap<float>> phoChargedWorstVtxIsoToken_;
0160 edm::EDGetTokenT<edm::ValueMap<float>> phoChargedWorstVtxGeomVetoIsoToken_;
0161 edm::EDGetTokenT<edm::ValueMap<float>> phoChargedPFPVIsoToken_;
0162
0163 edm::EDGetTokenT<edm::ValueMap<float>> phoPFECALClusIsolationToken_;
0164 edm::EDGetTokenT<edm::ValueMap<float>> phoPFHCALClusIsolationToken_;
0165
0166 const EcalClusterLazyTools::ESGetTokens ecalClusterESGetTokens_;
0167
0168 std::string valueMapPFCandPhoton_;
0169
0170 std::unique_ptr<PhotonIsolationCalculator> photonIsoCalculator_ = nullptr;
0171
0172
0173
0174
0175 std::vector<int> flagsexclEB_;
0176 std::vector<int> flagsexclEE_;
0177 std::vector<int> severitiesexclEB_;
0178 std::vector<int> severitiesexclEE_;
0179
0180 double multThresEB_;
0181 double multThresEE_;
0182 double hOverEConeSize_;
0183 bool checkHcalStatus_;
0184 double highEt_;
0185 double minR9Barrel_;
0186 double minR9Endcap_;
0187 bool runMIPTagger_;
0188 bool runMVABasedHaloTagger_;
0189
0190 RecoStepInfo recoStep_;
0191
0192 bool usePrimaryVertex_;
0193
0194 CaloGeometry const* caloGeom_ = nullptr;
0195
0196
0197 std::unique_ptr<PhotonMIPHaloTagger> photonMIPHaloTagger_ = nullptr;
0198
0199 std::unique_ptr<PhotonMVABasedHaloTagger> photonMVABasedHaloTagger_ = nullptr;
0200
0201 std::vector<double> preselCutValuesBarrel_;
0202 std::vector<double> preselCutValuesEndcap_;
0203
0204 std::unique_ptr<PhotonEnergyCorrector> photonEnergyCorrector_ = nullptr;
0205 std::string candidateP4type_;
0206
0207 const edm::ESGetToken<CaloTopology, CaloTopologyRecord> caloTopologyToken_;
0208 const edm::ESGetToken<CaloGeometry, CaloGeometryRecord> caloGeometryToken_;
0209 const edm::ESGetToken<EcalPFRecHitThresholds, EcalPFRecHitThresholdsRcd> ecalPFRechitThresholdsToken_;
0210
0211
0212 std::unique_ptr<ElectronHcalHelper> hcalHelperCone_;
0213 std::unique_ptr<ElectronHcalHelper> hcalHelperBc_;
0214 bool hcalRun2EffDepth_;
0215
0216
0217 bool dnnPFidEnabled_;
0218 std::vector<tensorflow::Session*> tfSessions_;
0219
0220 double ecaldrMax_;
0221 double ecaldrVetoBarrel_;
0222 double ecaldrVetoEndcap_;
0223 double ecaletaStripBarrel_;
0224 double ecaletaStripEndcap_;
0225 double ecalenergyBarrel_;
0226 double ecalenergyEndcap_;
0227 typedef EcalPFClusterIsolation<reco::Photon> PhotonEcalPFClusterIsolation;
0228 std::unique_ptr<PhotonEcalPFClusterIsolation> ecalisoAlgo = nullptr;
0229 edm::EDGetTokenT<reco::PFClusterCollection> pfClusterProducer_;
0230
0231 bool useHF_;
0232 double hcaldrMax_;
0233 double hcaldrVetoBarrel_;
0234 double hcaldrVetoEndcap_;
0235 double hcaletaStripBarrel_;
0236 double hcaletaStripEndcap_;
0237 double hcalenergyBarrel_;
0238 double hcalenergyEndcap_;
0239 double hcaluseEt_;
0240
0241 edm::EDGetTokenT<reco::PFClusterCollection> pfClusterProducerHCAL_;
0242 edm::EDGetTokenT<reco::PFClusterCollection> pfClusterProducerHFEM_;
0243 edm::EDGetTokenT<reco::PFClusterCollection> pfClusterProducerHFHAD_;
0244
0245 typedef HcalPFClusterIsolation<reco::Photon> PhotonHcalPFClusterIsolation;
0246 std::unique_ptr<PhotonHcalPFClusterIsolation> hcalisoAlgo = nullptr;
0247 };
0248
0249 #include "FWCore/Framework/interface/MakerMacros.h"
0250 DEFINE_FWK_MODULE(GEDPhotonProducer);
0251
0252 namespace {
0253 inline double ptFast(const double energy, const math::XYZPoint& position, const math::XYZPoint& origin) {
0254 const auto v = position - origin;
0255 return energy * std::sqrt(v.perp2() / v.mag2());
0256 }
0257 }
0258
0259 GEDPhotonProducer::RecoStepInfo::RecoStepInfo(const std::string& step) : flags_(0) {
0260 if (step == "final")
0261 flags_ = kFinal;
0262 else if (step == "oot")
0263 flags_ = kOOT;
0264 else if (step == "ootfinal")
0265 flags_ = (kOOT | kFinal);
0266 else if (step == "tmp")
0267 flags_ = 0;
0268 else {
0269 throw cms::Exception("InvalidConfig")
0270 << " reconstructStep " << step << " is invalid, the options are: tmp, final,oot or ootfinal" << std::endl;
0271 }
0272 }
0273
0274 GEDPhotonProducer::GEDPhotonProducer(const edm::ParameterSet& config, const CacheData* gcache)
0275 : photonProducer_{config.getParameter<edm::InputTag>("photonProducer")},
0276 ecalClusterESGetTokens_{consumesCollector()},
0277 recoStep_(config.getParameter<std::string>("reconstructionStep")),
0278 caloTopologyToken_{esConsumes()},
0279 caloGeometryToken_{esConsumes()},
0280 ecalPFRechitThresholdsToken_{esConsumes()},
0281 hcalHelperCone_(nullptr),
0282 hcalHelperBc_(nullptr) {
0283 if (recoStep_.isFinal()) {
0284 photonProducerT_ = consumes(photonProducer_);
0285 pfCandidates_ = consumes(config.getParameter<edm::InputTag>("pfCandidates"));
0286
0287 const edm::ParameterSet& pfIsolCfg = config.getParameter<edm::ParameterSet>("pfIsolCfg");
0288 auto getVMToken = [&pfIsolCfg, this](const std::string& name) {
0289 return consumes(pfIsolCfg.getParameter<edm::InputTag>(name));
0290 };
0291 phoChargedIsolationToken_ = getVMToken("chargedHadronIso");
0292 phoNeutralHadronIsolationToken_ = getVMToken("neutralHadronIso");
0293 phoPhotonIsolationToken_ = getVMToken("photonIso");
0294 phoChargedWorstVtxIsoToken_ = getVMToken("chargedHadronWorstVtxIso");
0295 phoChargedWorstVtxGeomVetoIsoToken_ = getVMToken("chargedHadronWorstVtxGeomVetoIso");
0296 phoChargedPFPVIsoToken_ = getVMToken("chargedHadronPFPVIso");
0297
0298
0299
0300 if (config.exists("pfECALClusIsolation")) {
0301 phoPFECALClusIsolationToken_ = consumes(config.getParameter<edm::InputTag>("pfECALClusIsolation"));
0302 }
0303 if (config.exists("pfHCALClusIsolation")) {
0304 phoPFHCALClusIsolationToken_ = consumes(config.getParameter<edm::InputTag>("pfHCALClusIsolation"));
0305 }
0306
0307 } else {
0308 photonCoreProducerT_ = consumes(photonProducer_);
0309 }
0310
0311 auto pfEg = config.getParameter<edm::InputTag>("pfEgammaCandidates");
0312 if (not pfEg.label().empty()) {
0313 pfEgammaCandidates_ = consumes(pfEg);
0314 }
0315 barrelEcalHits_ = consumes(config.getParameter<edm::InputTag>("barrelEcalHits"));
0316 endcapEcalHits_ = consumes(config.getParameter<edm::InputTag>("endcapEcalHits"));
0317 preshowerHits_ = consumes(config.getParameter<edm::InputTag>("preshowerHits"));
0318 vertexProducer_ = consumes(config.getParameter<edm::InputTag>("primaryVertexProducer"));
0319
0320 auto hbhetag = config.getParameter<edm::InputTag>("hbheRecHits");
0321 if (not hbhetag.label().empty())
0322 hbheRecHits_ = consumes<HBHERecHitCollection>(hbhetag);
0323
0324
0325 photonCollection_ = config.getParameter<std::string>("outputPhotonCollection");
0326 multThresEB_ = config.getParameter<double>("multThresEB");
0327 multThresEE_ = config.getParameter<double>("multThresEE");
0328 hOverEConeSize_ = config.getParameter<double>("hOverEConeSize");
0329 highEt_ = config.getParameter<double>("highEt");
0330
0331 minR9Barrel_ = config.getParameter<double>("minR9Barrel");
0332 minR9Endcap_ = config.getParameter<double>("minR9Endcap");
0333 usePrimaryVertex_ = config.getParameter<bool>("usePrimaryVertex");
0334 runMIPTagger_ = config.getParameter<bool>("runMIPTagger");
0335 runMVABasedHaloTagger_ = config.getParameter<bool>("runMVABasedHaloTagger");
0336
0337 candidateP4type_ = config.getParameter<std::string>("candidateP4type");
0338 valueMapPFCandPhoton_ = config.getParameter<std::string>("valueMapPhotons");
0339
0340
0341
0342 auto const& flagnamesEB = config.getParameter<std::vector<std::string>>("RecHitFlagToBeExcludedEB");
0343 auto const& flagnamesEE = config.getParameter<std::vector<std::string>>("RecHitFlagToBeExcludedEE");
0344
0345 flagsexclEB_ = StringToEnumValue<EcalRecHit::Flags>(flagnamesEB);
0346 flagsexclEE_ = StringToEnumValue<EcalRecHit::Flags>(flagnamesEE);
0347
0348 auto const& severitynamesEB = config.getParameter<std::vector<std::string>>("RecHitSeverityToBeExcludedEB");
0349 auto const& severitynamesEE = config.getParameter<std::vector<std::string>>("RecHitSeverityToBeExcludedEE");
0350
0351 severitiesexclEB_ = StringToEnumValue<EcalSeverityLevel::SeverityLevel>(severitynamesEB);
0352 severitiesexclEE_ = StringToEnumValue<EcalSeverityLevel::SeverityLevel>(severitynamesEE);
0353
0354 photonEnergyCorrector_ = std::make_unique<PhotonEnergyCorrector>(config, consumesCollector());
0355
0356 checkHcalStatus_ = config.getParameter<bool>("checkHcalStatus");
0357 if (not hbheRecHits_.isUninitialized()) {
0358 ElectronHcalHelper::Configuration cfgCone, cfgBc;
0359 cfgCone.hOverEConeSize = hOverEConeSize_;
0360 if (cfgCone.hOverEConeSize > 0) {
0361 cfgCone.onlyBehindCluster = false;
0362 cfgCone.checkHcalStatus = checkHcalStatus_;
0363
0364 cfgCone.hbheRecHits = hbheRecHits_;
0365
0366 cfgCone.eThresHB = config.getParameter<EgammaHcalIsolation::arrayHB>("recHitEThresholdHB");
0367 cfgCone.maxSeverityHB = config.getParameter<int>("maxHcalRecHitSeverity");
0368 cfgCone.eThresHE = config.getParameter<EgammaHcalIsolation::arrayHE>("recHitEThresholdHE");
0369 cfgCone.maxSeverityHE = cfgCone.maxSeverityHB;
0370 }
0371 cfgBc.hOverEConeSize = 0.;
0372 cfgBc.onlyBehindCluster = true;
0373 cfgBc.checkHcalStatus = checkHcalStatus_;
0374
0375 cfgBc.hbheRecHits = hbheRecHits_;
0376
0377 cfgBc.eThresHB = config.getParameter<EgammaHcalIsolation::arrayHB>("recHitEThresholdHB");
0378 cfgBc.maxSeverityHB = config.getParameter<int>("maxHcalRecHitSeverity");
0379 cfgBc.eThresHE = config.getParameter<EgammaHcalIsolation::arrayHE>("recHitEThresholdHE");
0380 cfgBc.maxSeverityHE = cfgBc.maxSeverityHB;
0381
0382 hcalHelperCone_ = std::make_unique<ElectronHcalHelper>(cfgCone, consumesCollector());
0383 hcalHelperBc_ = std::make_unique<ElectronHcalHelper>(cfgBc, consumesCollector());
0384 }
0385
0386 hcalRun2EffDepth_ = config.getParameter<bool>("hcalRun2EffDepth");
0387
0388
0389 preselCutValuesBarrel_ = {config.getParameter<double>("minSCEtBarrel"),
0390 config.getParameter<double>("maxHoverEBarrel"),
0391 config.getParameter<double>("ecalRecHitSumEtOffsetBarrel"),
0392 config.getParameter<double>("ecalRecHitSumEtSlopeBarrel"),
0393 config.getParameter<double>("hcalRecHitSumEtOffsetBarrel"),
0394 config.getParameter<double>("hcalRecHitSumEtSlopeBarrel"),
0395 config.getParameter<double>("nTrackSolidConeBarrel"),
0396 config.getParameter<double>("nTrackHollowConeBarrel"),
0397 config.getParameter<double>("trackPtSumSolidConeBarrel"),
0398 config.getParameter<double>("trackPtSumHollowConeBarrel"),
0399 config.getParameter<double>("sigmaIetaIetaCutBarrel")};
0400
0401 preselCutValuesEndcap_ = {config.getParameter<double>("minSCEtEndcap"),
0402 config.getParameter<double>("maxHoverEEndcap"),
0403 config.getParameter<double>("ecalRecHitSumEtOffsetEndcap"),
0404 config.getParameter<double>("ecalRecHitSumEtSlopeEndcap"),
0405 config.getParameter<double>("hcalRecHitSumEtOffsetEndcap"),
0406 config.getParameter<double>("hcalRecHitSumEtSlopeEndcap"),
0407 config.getParameter<double>("nTrackSolidConeEndcap"),
0408 config.getParameter<double>("nTrackHollowConeEndcap"),
0409 config.getParameter<double>("trackPtSumSolidConeEndcap"),
0410 config.getParameter<double>("trackPtSumHollowConeEndcap"),
0411 config.getParameter<double>("sigmaIetaIetaCutEndcap")};
0412
0413
0414
0415 if (!recoStep_.isFinal()) {
0416 photonIsoCalculator_ = std::make_unique<PhotonIsolationCalculator>();
0417 edm::ParameterSet isolationSumsCalculatorSet = config.getParameter<edm::ParameterSet>("isolationSumsCalculatorSet");
0418 photonIsoCalculator_->setup(isolationSumsCalculatorSet,
0419 flagsexclEB_,
0420 flagsexclEE_,
0421 severitiesexclEB_,
0422 severitiesexclEE_,
0423 consumesCollector());
0424 photonMIPHaloTagger_ = std::make_unique<PhotonMIPHaloTagger>();
0425 edm::ParameterSet mipVariableSet = config.getParameter<edm::ParameterSet>("mipVariableSet");
0426 photonMIPHaloTagger_->setup(mipVariableSet, consumesCollector());
0427 }
0428
0429 if (recoStep_.isFinal() && runMVABasedHaloTagger_) {
0430 edm::ParameterSet mvaBasedHaloVariableSet = config.getParameter<edm::ParameterSet>("mvaBasedHaloVariableSet");
0431 photonMVABasedHaloTagger_ =
0432 std::make_unique<PhotonMVABasedHaloTagger>(mvaBasedHaloVariableSet, consumesCollector());
0433 }
0434
0435
0436 const edm::ParameterSet& pfECALClusIsolCfg = config.getParameter<edm::ParameterSet>("pfECALClusIsolCfg");
0437 pfClusterProducer_ =
0438 consumes<reco::PFClusterCollection>(pfECALClusIsolCfg.getParameter<edm::InputTag>("pfClusterProducer"));
0439 ecaldrMax_ = pfECALClusIsolCfg.getParameter<double>("drMax");
0440 ecaldrVetoBarrel_ = pfECALClusIsolCfg.getParameter<double>("drVetoBarrel");
0441 ecaldrVetoEndcap_ = pfECALClusIsolCfg.getParameter<double>("drVetoEndcap");
0442 ecaletaStripBarrel_ = pfECALClusIsolCfg.getParameter<double>("etaStripBarrel");
0443 ecaletaStripEndcap_ = pfECALClusIsolCfg.getParameter<double>("etaStripEndcap");
0444 ecalenergyBarrel_ = pfECALClusIsolCfg.getParameter<double>("energyBarrel");
0445 ecalenergyEndcap_ = pfECALClusIsolCfg.getParameter<double>("energyEndcap");
0446
0447 const edm::ParameterSet& pfHCALClusIsolCfg = config.getParameter<edm::ParameterSet>("pfHCALClusIsolCfg");
0448 pfClusterProducerHCAL_ = consumes(pfHCALClusIsolCfg.getParameter<edm::InputTag>("pfClusterProducerHCAL"));
0449 pfClusterProducerHFEM_ = consumes(pfHCALClusIsolCfg.getParameter<edm::InputTag>("pfClusterProducerHFEM"));
0450 pfClusterProducerHFHAD_ = consumes(pfHCALClusIsolCfg.getParameter<edm::InputTag>("pfClusterProducerHFHAD"));
0451 useHF_ = pfHCALClusIsolCfg.getParameter<bool>("useHF");
0452 hcaldrMax_ = pfHCALClusIsolCfg.getParameter<double>("drMax");
0453 hcaldrVetoBarrel_ = pfHCALClusIsolCfg.getParameter<double>("drVetoBarrel");
0454 hcaldrVetoEndcap_ = pfHCALClusIsolCfg.getParameter<double>("drVetoEndcap");
0455 hcaletaStripBarrel_ = pfHCALClusIsolCfg.getParameter<double>("etaStripBarrel");
0456 hcaletaStripEndcap_ = pfHCALClusIsolCfg.getParameter<double>("etaStripEndcap");
0457 hcalenergyBarrel_ = pfHCALClusIsolCfg.getParameter<double>("energyBarrel");
0458 hcalenergyEndcap_ = pfHCALClusIsolCfg.getParameter<double>("energyEndcap");
0459 hcaluseEt_ = pfHCALClusIsolCfg.getParameter<bool>("useEt");
0460
0461
0462 produces<reco::PhotonCollection>(photonCollection_);
0463 if (not pfEgammaCandidates_.isUninitialized()) {
0464 produces<edm::ValueMap<reco::PhotonRef>>(valueMapPFCandPhoton_);
0465 }
0466
0467 const auto& pset_dnn = config.getParameter<edm::ParameterSet>("PhotonDNNPFid");
0468 dnnPFidEnabled_ = pset_dnn.getParameter<bool>("enabled");
0469 if (dnnPFidEnabled_) {
0470 tfSessions_ = gcache->photonDNNEstimator->getSessions();
0471 }
0472 }
0473
0474 std::unique_ptr<CacheData> GEDPhotonProducer::initializeGlobalCache(const edm::ParameterSet& config) {
0475
0476 return std::make_unique<CacheData>(config);
0477 }
0478
0479 void GEDPhotonProducer::endStream() {
0480 for (auto session : tfSessions_) {
0481 tensorflow::closeSession(session);
0482 }
0483 }
0484
0485 void GEDPhotonProducer::produce(edm::Event& theEvent, const edm::EventSetup& eventSetup) {
0486 using namespace edm;
0487
0488 auto outputPhotonCollection_p = std::make_unique<reco::PhotonCollection>();
0489 edm::ValueMap<reco::PhotonRef> pfEGCandToPhotonMap;
0490
0491
0492 bool validPhotonCoreHandle = false;
0493 Handle<reco::PhotonCoreCollection> photonCoreHandle;
0494 bool validPhotonHandle = false;
0495 Handle<reco::PhotonCollection> photonHandle;
0496
0497 edm::Handle<edm::ValueMap<float>> phoChargedIsolationMap;
0498 edm::Handle<edm::ValueMap<float>> phoNeutralHadronIsolationMap;
0499 edm::Handle<edm::ValueMap<float>> phoPhotonIsolationMap;
0500 edm::Handle<edm::ValueMap<float>> phoChargedWorstVtxIsoMap;
0501 edm::Handle<edm::ValueMap<float>> phoChargedWorstVtxGeomVetoIsoMap;
0502 edm::Handle<edm::ValueMap<float>> phoChargedPFPVIsoMap;
0503
0504 edm::Handle<edm::ValueMap<float>> phoPFECALClusIsolationMap;
0505 edm::Handle<edm::ValueMap<float>> phoPFHCALClusIsolationMap;
0506
0507 if (recoStep_.isFinal()) {
0508 theEvent.getByToken(photonProducerT_, photonHandle);
0509
0510 theEvent.getByToken(phoChargedIsolationToken_, phoChargedIsolationMap);
0511 theEvent.getByToken(phoNeutralHadronIsolationToken_, phoNeutralHadronIsolationMap);
0512 theEvent.getByToken(phoPhotonIsolationToken_, phoPhotonIsolationMap);
0513 theEvent.getByToken(phoChargedWorstVtxIsoToken_, phoChargedWorstVtxIsoMap);
0514 theEvent.getByToken(phoChargedWorstVtxGeomVetoIsoToken_, phoChargedWorstVtxGeomVetoIsoMap);
0515 theEvent.getByToken(phoChargedPFPVIsoToken_, phoChargedPFPVIsoMap);
0516
0517
0518 if (!phoPFECALClusIsolationToken_.isUninitialized()) {
0519 theEvent.getByToken(phoPFECALClusIsolationToken_, phoPFECALClusIsolationMap);
0520 }
0521 if (!phoPFHCALClusIsolationToken_.isUninitialized()) {
0522 theEvent.getByToken(phoPFHCALClusIsolationToken_, phoPFHCALClusIsolationMap);
0523 }
0524
0525 if (photonHandle.isValid()) {
0526 validPhotonHandle = true;
0527 } else {
0528 throw cms::Exception("GEDPhotonProducer") << "Error! Can't get the product " << photonProducer_.label() << "\n";
0529 }
0530 } else {
0531 theEvent.getByToken(photonCoreProducerT_, photonCoreHandle);
0532 if (photonCoreHandle.isValid()) {
0533 validPhotonCoreHandle = true;
0534 } else {
0535 throw cms::Exception("GEDPhotonProducer")
0536 << "Error! Can't get the photonCoreProducer " << photonProducer_.label() << "\n";
0537 }
0538 }
0539
0540
0541 auto const& barrelRecHits = theEvent.get(barrelEcalHits_);
0542 auto const& endcapRecHits = theEvent.get(endcapEcalHits_);
0543 auto const& preshowerRecHits = theEvent.get(preshowerHits_);
0544
0545 Handle<reco::PFCandidateCollection> pfEGCandidateHandle;
0546
0547 if (not pfEgammaCandidates_.isUninitialized()) {
0548 theEvent.getByToken(pfEgammaCandidates_, pfEGCandidateHandle);
0549 if (!pfEGCandidateHandle.isValid()) {
0550 throw cms::Exception("GEDPhotonProducer") << "Error! Can't get the pfEgammaCandidates";
0551 }
0552 }
0553
0554 Handle<reco::PFCandidateCollection> pfCandidateHandle;
0555
0556 if (recoStep_.isFinal()) {
0557
0558 theEvent.getByToken(pfCandidates_, pfCandidateHandle);
0559
0560 if (!pfCandidateHandle.isValid() && !recoStep_.isOOT()) {
0561 throw cms::Exception("GEDPhotonProducer") << "Error! Can't get the pfCandidates";
0562 }
0563 }
0564
0565
0566 caloGeom_ = &eventSetup.getData(caloGeometryToken_);
0567
0568
0569 if (hcalHelperCone_ != nullptr and hcalHelperBc_ != nullptr) {
0570 hcalHelperCone_->beginEvent(theEvent, eventSetup);
0571 hcalHelperBc_->beginEvent(theEvent, eventSetup);
0572 }
0573
0574 auto const& topology = eventSetup.getData(caloTopologyToken_);
0575 auto const& thresholds = eventSetup.getData(ecalPFRechitThresholdsToken_);
0576
0577
0578 const reco::VertexCollection dummyVC;
0579 auto const& vertexCollection{usePrimaryVertex_ ? theEvent.get(vertexProducer_) : dummyVC};
0580
0581
0582
0583
0584
0585 photonEnergyCorrector_->init(eventSetup);
0586 if (photonEnergyCorrector_->gedRegression()) {
0587 photonEnergyCorrector_->gedRegression()->setEvent(theEvent);
0588 photonEnergyCorrector_->gedRegression()->setEventContent(eventSetup);
0589 }
0590
0591
0592 ecalisoAlgo = std::make_unique<PhotonEcalPFClusterIsolation>(ecaldrMax_,
0593 ecaldrVetoBarrel_,
0594 ecaldrVetoEndcap_,
0595 ecaletaStripBarrel_,
0596 ecaletaStripEndcap_,
0597 ecalenergyBarrel_,
0598 ecalenergyEndcap_);
0599
0600 hcalisoAlgo = std::make_unique<PhotonHcalPFClusterIsolation>(hcaldrMax_,
0601 hcaldrVetoBarrel_,
0602 hcaldrVetoEndcap_,
0603 hcaletaStripBarrel_,
0604 hcaletaStripEndcap_,
0605 hcalenergyBarrel_,
0606 hcalenergyEndcap_,
0607 hcaluseEt_);
0608
0609 int iSC = 0;
0610
0611 if (validPhotonCoreHandle)
0612 fillPhotonCollection(theEvent,
0613 eventSetup,
0614 photonCoreHandle,
0615 &topology,
0616 &barrelRecHits,
0617 &endcapRecHits,
0618 &preshowerRecHits,
0619 hcalHelperCone_.get(),
0620 hcalHelperBc_.get(),
0621
0622 vertexCollection,
0623 *outputPhotonCollection_p,
0624 iSC,
0625 thresholds);
0626
0627 iSC = 0;
0628 if (validPhotonHandle && recoStep_.isFinal())
0629 fillPhotonCollection(theEvent,
0630 eventSetup,
0631 photonHandle,
0632 pfCandidateHandle,
0633 pfEGCandidateHandle,
0634 theEvent.get(vertexProducer_),
0635 *outputPhotonCollection_p,
0636 iSC,
0637 phoChargedIsolationMap,
0638 phoNeutralHadronIsolationMap,
0639 phoPhotonIsolationMap,
0640 phoChargedWorstVtxIsoMap,
0641 phoChargedWorstVtxGeomVetoIsoMap,
0642 phoChargedPFPVIsoMap,
0643 phoPFECALClusIsolationMap,
0644 phoPFHCALClusIsolationMap);
0645
0646
0647 edm::LogInfo("GEDPhotonProducer") << " Put in the event " << iSC << " Photon Candidates \n";
0648
0649
0650 if (hcalRun2EffDepth_) {
0651 for (auto& pho : *outputPhotonCollection_p)
0652 pho.hcalToRun2EffDepth();
0653 }
0654 const auto photonOrphHandle = theEvent.put(std::move(outputPhotonCollection_p), photonCollection_);
0655
0656 if (!recoStep_.isFinal() && not pfEgammaCandidates_.isUninitialized()) {
0657
0658 auto pfEGCandToPhotonMap_p = std::make_unique<edm::ValueMap<reco::PhotonRef>>();
0659 edm::ValueMap<reco::PhotonRef>::Filler filler(*pfEGCandToPhotonMap_p);
0660 unsigned nObj = pfEGCandidateHandle->size();
0661 std::vector<reco::PhotonRef> values(nObj);
0662
0663 for (unsigned int lCand = 0; lCand < nObj; lCand++) {
0664 reco::PFCandidateRef pfCandRef(reco::PFCandidateRef(pfEGCandidateHandle, lCand));
0665 reco::SuperClusterRef pfScRef = pfCandRef->superClusterRef();
0666
0667 for (unsigned int lSC = 0; lSC < photonOrphHandle->size(); lSC++) {
0668 reco::PhotonRef photonRef(reco::PhotonRef(photonOrphHandle, lSC));
0669 reco::SuperClusterRef scRef = photonRef->superCluster();
0670 if (pfScRef != scRef)
0671 continue;
0672 values[lCand] = photonRef;
0673 }
0674 }
0675
0676 filler.insert(pfEGCandidateHandle, values.begin(), values.end());
0677 filler.fill();
0678 theEvent.put(std::move(pfEGCandToPhotonMap_p), valueMapPFCandPhoton_);
0679 }
0680 }
0681
0682 void GEDPhotonProducer::fillPhotonCollection(edm::Event& evt,
0683 edm::EventSetup const& es,
0684 const edm::Handle<reco::PhotonCoreCollection>& photonCoreHandle,
0685 const CaloTopology* topology,
0686 const EcalRecHitCollection* ecalBarrelHits,
0687 const EcalRecHitCollection* ecalEndcapHits,
0688 const EcalRecHitCollection* preshowerHits,
0689 const ElectronHcalHelper* hcalHelperCone,
0690 const ElectronHcalHelper* hcalHelperBc,
0691 const reco::VertexCollection& vertexCollection,
0692 reco::PhotonCollection& outputPhotonCollection,
0693 int& iSC,
0694 EcalPFRecHitThresholds const& thresholds) {
0695 const EcalRecHitCollection* hits = nullptr;
0696 std::vector<double> preselCutValues;
0697 std::vector<int> flags_, severitiesexcl_;
0698
0699 for (unsigned int lSC = 0; lSC < photonCoreHandle->size(); lSC++) {
0700 reco::PhotonCoreRef coreRef(reco::PhotonCoreRef(photonCoreHandle, lSC));
0701 reco::SuperClusterRef parentSCRef = coreRef->parentSuperCluster();
0702 reco::SuperClusterRef scRef = coreRef->superCluster();
0703
0704
0705 iSC++;
0706
0707 DetId::Detector thedet = scRef->seed()->hitsAndFractions()[0].first.det();
0708 int subdet = scRef->seed()->hitsAndFractions()[0].first.subdetId();
0709 if (subdet == EcalBarrel) {
0710 preselCutValues = preselCutValuesBarrel_;
0711 hits = ecalBarrelHits;
0712 flags_ = flagsexclEB_;
0713 severitiesexcl_ = severitiesexclEB_;
0714 } else if (subdet == EcalEndcap) {
0715 preselCutValues = preselCutValuesEndcap_;
0716 hits = ecalEndcapHits;
0717 flags_ = flagsexclEE_;
0718 severitiesexcl_ = severitiesexclEE_;
0719 } else if (EcalTools::isHGCalDet(thedet)) {
0720 preselCutValues = preselCutValuesEndcap_;
0721 hits = nullptr;
0722 flags_ = flagsexclEE_;
0723 severitiesexcl_ = severitiesexclEE_;
0724 } else {
0725 edm::LogWarning("") << "GEDPhotonProducer: do not know if it is a barrel or endcap SuperCluster: " << thedet
0726 << ' ' << subdet;
0727 }
0728
0729
0730 if (parentSCRef.isNonnull() &&
0731 ptFast(parentSCRef->energy(), parentSCRef->position(), {0, 0, 0}) <= preselCutValues[0])
0732 continue;
0733
0734 float maxXtal = (hits != nullptr ? EcalClusterTools::eMax(*(scRef->seed()), hits) : 0.f);
0735
0736
0737
0738 float e1x5 = (hits != nullptr ? EcalClusterTools::e1x5(*(scRef->seed()), hits, topology) : 0.f);
0739 float e2x5 = (hits != nullptr ? EcalClusterTools::e2x5Max(*(scRef->seed()), hits, topology) : 0.f);
0740 float e3x3 = (hits != nullptr ? EcalClusterTools::e3x3(*(scRef->seed()), hits, topology) : 0.f);
0741 float e5x5 = (hits != nullptr ? EcalClusterTools::e5x5(*(scRef->seed()), hits, topology) : 0.f);
0742 const auto& cov = (hits != nullptr ? EcalClusterTools::covariances(*(scRef->seed()), hits, topology, caloGeom_)
0743 : std::array<float, 3>({{0.f, 0.f, 0.f}}));
0744
0745 const auto& locCov = (hits != nullptr ? EcalClusterTools::localCovariances(*(scRef->seed()), hits, topology)
0746 : std::array<float, 3>({{0.f, 0.f, 0.f}}));
0747
0748 float sigmaEtaEta = std::sqrt(cov[0]);
0749 float sigmaIetaIeta = std::sqrt(locCov[0]);
0750
0751 float full5x5_maxXtal = (hits != nullptr ? noZS::EcalClusterTools::eMax(*(scRef->seed()), hits) : 0.f);
0752
0753
0754 float full5x5_e1x5 = (hits != nullptr ? noZS::EcalClusterTools::e1x5(*(scRef->seed()), hits, topology) : 0.f);
0755 float full5x5_e2x5 = (hits != nullptr ? noZS::EcalClusterTools::e2x5Max(*(scRef->seed()), hits, topology) : 0.f);
0756 float full5x5_e3x3 = (hits != nullptr ? noZS::EcalClusterTools::e3x3(*(scRef->seed()), hits, topology) : 0.f);
0757 float full5x5_e5x5 = (hits != nullptr ? noZS::EcalClusterTools::e5x5(*(scRef->seed()), hits, topology) : 0.f);
0758 const auto& full5x5_cov =
0759 (hits != nullptr ? noZS::EcalClusterTools::covariances(*(scRef->seed()), hits, topology, caloGeom_)
0760 : std::array<float, 3>({{0.f, 0.f, 0.f}}));
0761
0762
0763
0764 const auto& full5x5_locCov =
0765 (hits != nullptr ? noZS::EcalClusterTools::localCovariances(*(scRef->seed()),
0766 hits,
0767 topology,
0768 EgammaLocalCovParamDefaults::kRelEnCut,
0769 &thresholds,
0770 multThresEB_,
0771 multThresEE_)
0772 : std::array<float, 3>({{0.f, 0.f, 0.f}}));
0773
0774 float full5x5_sigmaEtaEta = sqrt(full5x5_cov[0]);
0775 float full5x5_sigmaIetaIeta = sqrt(full5x5_locCov[0]);
0776 float full5x5_sigmaIetaIphi = full5x5_locCov[1];
0777
0778
0779 math::XYZPoint caloPosition = scRef->position();
0780
0781
0782 double photonEnergy = 1.;
0783 math::XYZPoint vtx(0., 0., 0.);
0784 if (!vertexCollection.empty())
0785 vtx = vertexCollection.begin()->position();
0786
0787 math::XYZVector direction = caloPosition - vtx;
0788
0789 math::XYZVector momentum = direction.unit();
0790
0791
0792 math::XYZTLorentzVectorD p4(momentum.x(), momentum.y(), momentum.z(), photonEnergy);
0793 reco::Photon newCandidate(p4, caloPosition, coreRef, vtx);
0794
0795
0796
0797
0798
0799 reco::Photon::FiducialFlags fiducialFlags;
0800 reco::Photon::IsolationVariables isolVarR03, isolVarR04;
0801 if (!EcalTools::isHGCalDet(thedet)) {
0802 photonIsoCalculator_->calculate(&newCandidate, evt, es, fiducialFlags, isolVarR04, isolVarR03);
0803 }
0804 newCandidate.setFiducialVolumeFlags(fiducialFlags);
0805 newCandidate.setIsolationVariables(isolVarR04, isolVarR03);
0806
0807
0808 reco::Photon::ShowerShape showerShape;
0809 showerShape.e1x5 = e1x5;
0810 showerShape.e2x5 = e2x5;
0811 showerShape.e3x3 = e3x3;
0812 showerShape.e5x5 = e5x5;
0813 showerShape.maxEnergyXtal = maxXtal;
0814 showerShape.sigmaEtaEta = sigmaEtaEta;
0815 showerShape.sigmaIetaIeta = sigmaIetaIeta;
0816 for (uint id = 0; id < showerShape.hcalOverEcal.size(); ++id) {
0817 showerShape.hcalOverEcal[id] =
0818 (hcalHelperCone != nullptr) ? hcalHelperCone->hcalESum(*scRef, id + 1) / scRef->energy() : 0.f;
0819
0820 showerShape.hcalOverEcalBc[id] =
0821 (hcalHelperBc != nullptr) ? hcalHelperBc->hcalESum(*scRef, id + 1) / scRef->energy() : 0.f;
0822 }
0823 showerShape.invalidHcal = (hcalHelperBc != nullptr) ? !hcalHelperBc->hasActiveHcal(*scRef) : false;
0824 if (hcalHelperBc != nullptr)
0825 showerShape.hcalTowersBehindClusters = hcalHelperBc->hcalTowersBehindClusters(*scRef);
0826 showerShape.pre7DepthHcal = false;
0827
0828
0829 const float spp = (!edm::isFinite(locCov[2]) ? 0. : sqrt(locCov[2]));
0830 const float sep = locCov[1];
0831 showerShape.sigmaIetaIphi = sep;
0832 showerShape.sigmaIphiIphi = spp;
0833 showerShape.e2nd = (hits != nullptr ? EcalClusterTools::e2nd(*(scRef->seed()), hits) : 0.f);
0834 showerShape.eTop = (hits != nullptr ? EcalClusterTools::eTop(*(scRef->seed()), hits, topology) : 0.f);
0835 showerShape.eLeft = (hits != nullptr ? EcalClusterTools::eLeft(*(scRef->seed()), hits, topology) : 0.f);
0836 showerShape.eRight = (hits != nullptr ? EcalClusterTools::eRight(*(scRef->seed()), hits, topology) : 0.f);
0837 showerShape.eBottom = (hits != nullptr ? EcalClusterTools::eBottom(*(scRef->seed()), hits, topology) : 0.f);
0838 showerShape.e1x3 = (hits != nullptr ? EcalClusterTools::e1x3(*(scRef->seed()), hits, topology) : 0.f);
0839 showerShape.e2x2 = (hits != nullptr ? EcalClusterTools::e2x2(*(scRef->seed()), hits, topology) : 0.f);
0840 showerShape.e2x5Max = (hits != nullptr ? EcalClusterTools::e2x5Max(*(scRef->seed()), hits, topology) : 0.f);
0841 showerShape.e2x5Left = (hits != nullptr ? EcalClusterTools::e2x5Left(*(scRef->seed()), hits, topology) : 0.f);
0842 showerShape.e2x5Right = (hits != nullptr ? EcalClusterTools::e2x5Right(*(scRef->seed()), hits, topology) : 0.f);
0843 showerShape.e2x5Top = (hits != nullptr ? EcalClusterTools::e2x5Top(*(scRef->seed()), hits, topology) : 0.f);
0844 showerShape.e2x5Bottom = (hits != nullptr ? EcalClusterTools::e2x5Bottom(*(scRef->seed()), hits, topology) : 0.f);
0845 if (hits) {
0846 Cluster2ndMoments clus2ndMoments = EcalClusterTools::cluster2ndMoments(*(scRef->seed()), *hits);
0847 showerShape.smMajor = clus2ndMoments.sMaj;
0848 showerShape.smMinor = clus2ndMoments.sMin;
0849 showerShape.smAlpha = clus2ndMoments.alpha;
0850 } else {
0851 showerShape.smMajor = 0.f;
0852 showerShape.smMinor = 0.f;
0853 showerShape.smAlpha = 0.f;
0854 }
0855
0856
0857 EcalClusterLazyTools toolsforES(
0858 evt, ecalClusterESGetTokens_.get(es), barrelEcalHits_, endcapEcalHits_, preshowerHits_);
0859 const float sigmaRR = toolsforES.eseffsirir(*scRef);
0860 showerShape.effSigmaRR = sigmaRR;
0861 newCandidate.setShowerShapeVariables(showerShape);
0862
0863 const reco::CaloCluster& seedCluster = *(scRef->seed());
0864 DetId seedXtalId = seedCluster.seed();
0865 int nSaturatedXtals = 0;
0866 bool isSeedSaturated = false;
0867 if (hits != nullptr) {
0868 const auto hitsAndFractions = scRef->hitsAndFractions();
0869 for (auto const& hitFractionPair : hitsAndFractions) {
0870 auto&& ecalRecHit = hits->find(hitFractionPair.first);
0871 if (ecalRecHit == hits->end())
0872 continue;
0873 if (ecalRecHit->checkFlag(EcalRecHit::Flags::kSaturated)) {
0874 nSaturatedXtals++;
0875 if (seedXtalId == ecalRecHit->detid())
0876 isSeedSaturated = true;
0877 }
0878 }
0879 }
0880 reco::Photon::SaturationInfo saturationInfo;
0881 saturationInfo.nSaturatedXtals = nSaturatedXtals;
0882 saturationInfo.isSeedSaturated = isSeedSaturated;
0883 newCandidate.setSaturationInfo(saturationInfo);
0884
0885
0886 reco::Photon::ShowerShape full5x5_showerShape;
0887 full5x5_showerShape.e1x5 = full5x5_e1x5;
0888 full5x5_showerShape.e2x5 = full5x5_e2x5;
0889 full5x5_showerShape.e3x3 = full5x5_e3x3;
0890 full5x5_showerShape.e5x5 = full5x5_e5x5;
0891 full5x5_showerShape.maxEnergyXtal = full5x5_maxXtal;
0892 full5x5_showerShape.sigmaEtaEta = full5x5_sigmaEtaEta;
0893 full5x5_showerShape.sigmaIetaIeta = full5x5_sigmaIetaIeta;
0894
0895 const float full5x5_spp = (!edm::isFinite(full5x5_locCov[2]) ? 0. : std::sqrt(full5x5_locCov[2]));
0896 const float full5x5_sep = full5x5_sigmaIetaIphi;
0897 full5x5_showerShape.sigmaIetaIphi = full5x5_sep;
0898 full5x5_showerShape.sigmaIphiIphi = full5x5_spp;
0899 full5x5_showerShape.e2nd = (hits != nullptr ? noZS::EcalClusterTools::e2nd(*(scRef->seed()), hits) : 0.f);
0900 full5x5_showerShape.eTop = (hits != nullptr ? noZS::EcalClusterTools::eTop(*(scRef->seed()), hits, topology) : 0.f);
0901 full5x5_showerShape.eLeft =
0902 (hits != nullptr ? noZS::EcalClusterTools::eLeft(*(scRef->seed()), hits, topology) : 0.f);
0903 full5x5_showerShape.eRight =
0904 (hits != nullptr ? noZS::EcalClusterTools::eRight(*(scRef->seed()), hits, topology) : 0.f);
0905 full5x5_showerShape.eBottom =
0906 (hits != nullptr ? noZS::EcalClusterTools::eBottom(*(scRef->seed()), hits, topology) : 0.f);
0907 full5x5_showerShape.e1x3 = (hits != nullptr ? noZS::EcalClusterTools::e1x3(*(scRef->seed()), hits, topology) : 0.f);
0908 full5x5_showerShape.e2x2 = (hits != nullptr ? noZS::EcalClusterTools::e2x2(*(scRef->seed()), hits, topology) : 0.f);
0909 full5x5_showerShape.e2x5Max =
0910 (hits != nullptr ? noZS::EcalClusterTools::e2x5Max(*(scRef->seed()), hits, topology) : 0.f);
0911 full5x5_showerShape.e2x5Left =
0912 (hits != nullptr ? noZS::EcalClusterTools::e2x5Left(*(scRef->seed()), hits, topology) : 0.f);
0913 full5x5_showerShape.e2x5Right =
0914 (hits != nullptr ? noZS::EcalClusterTools::e2x5Right(*(scRef->seed()), hits, topology) : 0.f);
0915 full5x5_showerShape.e2x5Top =
0916 (hits != nullptr ? noZS::EcalClusterTools::e2x5Top(*(scRef->seed()), hits, topology) : 0.f);
0917 full5x5_showerShape.e2x5Bottom =
0918 (hits != nullptr ? noZS::EcalClusterTools::e2x5Bottom(*(scRef->seed()), hits, topology) : 0.f);
0919 if (hits) {
0920 Cluster2ndMoments clus2ndMoments = noZS::EcalClusterTools::cluster2ndMoments(*(scRef->seed()), *hits);
0921 full5x5_showerShape.smMajor = clus2ndMoments.sMaj;
0922 full5x5_showerShape.smMinor = clus2ndMoments.sMin;
0923 full5x5_showerShape.smAlpha = clus2ndMoments.alpha;
0924 } else {
0925 full5x5_showerShape.smMajor = 0.f;
0926 full5x5_showerShape.smMinor = 0.f;
0927 full5x5_showerShape.smAlpha = 0.f;
0928 }
0929
0930 full5x5_showerShape.effSigmaRR = sigmaRR;
0931 for (uint id = 0; id < full5x5_showerShape.hcalOverEcal.size(); ++id) {
0932 full5x5_showerShape.hcalOverEcal[id] =
0933 (hcalHelperCone != nullptr) ? hcalHelperCone->hcalESum(*scRef, id + 1) / full5x5_e5x5 : 0.f;
0934 full5x5_showerShape.hcalOverEcalBc[id] =
0935 (hcalHelperBc != nullptr) ? hcalHelperBc->hcalESum(*scRef, id + 1) / full5x5_e5x5 : 0.f;
0936 }
0937 full5x5_showerShape.pre7DepthHcal = false;
0938 newCandidate.full5x5_setShowerShapeVariables(full5x5_showerShape);
0939
0940
0941 edm::Ptr<reco::PhotonCore> photonPtr(photonCoreHandle, lSC);
0942
0943
0944
0945
0946
0947 if (dnnPFidEnabled_) {
0948 auto clusterHandle = evt.getHandle(pfClusterProducer_);
0949 std::vector<edm::Handle<reco::PFClusterCollection>> clusterHandles{evt.getHandle(pfClusterProducerHCAL_)};
0950 if (useHF_) {
0951 clusterHandles.push_back(evt.getHandle(pfClusterProducerHFEM_));
0952 clusterHandles.push_back(evt.getHandle(pfClusterProducerHFHAD_));
0953 }
0954 reco::Photon::PflowIsolationVariables pfIso;
0955 pfIso.sumEcalClusterEt = ecalisoAlgo->getSum(newCandidate, clusterHandle);
0956 pfIso.sumHcalClusterEt = hcalisoAlgo->getSum(newCandidate, clusterHandles);
0957
0958 newCandidate.setPflowIsolationVariables(pfIso);
0959 }
0960
0961
0962
0963
0964
0965 if (!EcalTools::isHGCalDet(thedet)) {
0966 photonEnergyCorrector_->calculate(evt, newCandidate, subdet, vertexCollection, es);
0967 if (candidateP4type_ == "fromEcalEnergy") {
0968 newCandidate.setP4(newCandidate.p4(reco::Photon::ecal_photons));
0969 newCandidate.setCandidateP4type(reco::Photon::ecal_photons);
0970 newCandidate.setMass(0.0);
0971 } else if (candidateP4type_ == "fromRegression1") {
0972 newCandidate.setP4(newCandidate.p4(reco::Photon::regression1));
0973 newCandidate.setCandidateP4type(reco::Photon::regression1);
0974 newCandidate.setMass(0.0);
0975 } else if (candidateP4type_ == "fromRegression2") {
0976 newCandidate.setP4(newCandidate.p4(reco::Photon::regression2));
0977 newCandidate.setCandidateP4type(reco::Photon::regression2);
0978 newCandidate.setMass(0.0);
0979 } else if (candidateP4type_ == "fromRefinedSCRegression") {
0980 newCandidate.setP4(newCandidate.p4(reco::Photon::regression2));
0981 newCandidate.setCandidateP4type(reco::Photon::regression2);
0982 newCandidate.setMass(0.0);
0983 }
0984 } else {
0985 math::XYZVector gamma_momentum = direction.unit() * scRef->energy();
0986 math::PtEtaPhiMLorentzVector p4(gamma_momentum.rho(), gamma_momentum.eta(), gamma_momentum.phi(), 0.0);
0987 newCandidate.setP4(p4);
0988 newCandidate.setCandidateP4type(reco::Photon::ecal_photons);
0989
0990 reco::Photon::FiducialFlags fiducialFlags;
0991 fiducialFlags.isEE = true;
0992 newCandidate.setFiducialVolumeFlags(fiducialFlags);
0993 }
0994
0995
0996 reco::Photon::MIPVariables mipVar;
0997 if (subdet == EcalBarrel && runMIPTagger_) {
0998 photonMIPHaloTagger_->MIPcalculate(&newCandidate, evt, es, mipVar);
0999 newCandidate.setMIPVariables(mipVar);
1000 }
1001
1002
1003 bool isLooseEM = true;
1004 if (newCandidate.pt() < highEt_) {
1005 if (newCandidate.hadronicOverEm() >= preselCutValues[1])
1006 isLooseEM = false;
1007 if (newCandidate.ecalRecHitSumEtConeDR04() > preselCutValues[2] + preselCutValues[3] * newCandidate.pt())
1008 isLooseEM = false;
1009 if (newCandidate.hcalTowerSumEtConeDR04() > preselCutValues[4] + preselCutValues[5] * newCandidate.pt())
1010 isLooseEM = false;
1011 if (newCandidate.nTrkSolidConeDR04() > int(preselCutValues[6]))
1012 isLooseEM = false;
1013 if (newCandidate.nTrkHollowConeDR04() > int(preselCutValues[7]))
1014 isLooseEM = false;
1015 if (newCandidate.trkSumPtSolidConeDR04() > preselCutValues[8])
1016 isLooseEM = false;
1017 if (newCandidate.trkSumPtHollowConeDR04() > preselCutValues[9])
1018 isLooseEM = false;
1019 if (newCandidate.sigmaIetaIeta() > preselCutValues[10])
1020 isLooseEM = false;
1021 }
1022
1023 if (isLooseEM)
1024 outputPhotonCollection.push_back(newCandidate);
1025 }
1026
1027 if (dnnPFidEnabled_) {
1028
1029 LogDebug("GEDPhotonProducer") << "Getting DNN PFId for photons";
1030 const auto& dnn_photon_pfid = globalCache()->photonDNNEstimator->evaluate(outputPhotonCollection, tfSessions_);
1031 size_t ipho = 0;
1032 for (auto& photon : outputPhotonCollection) {
1033 const auto& [iModel, values] = dnn_photon_pfid[ipho];
1034 reco::Photon::PflowIDVariables pfID;
1035
1036 pfID.dnn = values[0];
1037 photon.setPflowIDVariables(pfID);
1038 ipho++;
1039 }
1040 }
1041 }
1042
1043 void GEDPhotonProducer::fillPhotonCollection(edm::Event& evt,
1044 edm::EventSetup const& es,
1045 const edm::Handle<reco::PhotonCollection>& photonHandle,
1046 const edm::Handle<reco::PFCandidateCollection> pfCandidateHandle,
1047 const edm::Handle<reco::PFCandidateCollection> pfEGCandidateHandle,
1048 reco::VertexCollection const& vertexCollection,
1049 reco::PhotonCollection& outputPhotonCollection,
1050 int& iSC,
1051 const edm::Handle<edm::ValueMap<float>>& chargedHadrons,
1052 const edm::Handle<edm::ValueMap<float>>& neutralHadrons,
1053 const edm::Handle<edm::ValueMap<float>>& photons,
1054 const edm::Handle<edm::ValueMap<float>>& chargedHadronsWorstVtx,
1055 const edm::Handle<edm::ValueMap<float>>& chargedHadronsWorstVtxGeomVeto,
1056 const edm::Handle<edm::ValueMap<float>>& chargedHadronsPFPV,
1057 const edm::Handle<edm::ValueMap<float>>& pfEcalClusters,
1058 const edm::Handle<edm::ValueMap<float>>& pfHcalClusters) {
1059 std::vector<double> preselCutValues;
1060
1061 for (unsigned int lSC = 0; lSC < photonHandle->size(); lSC++) {
1062 reco::PhotonRef phoRef(reco::PhotonRef(photonHandle, lSC));
1063 reco::SuperClusterRef parentSCRef = phoRef->parentSuperCluster();
1064 reco::SuperClusterRef scRef = phoRef->superCluster();
1065 DetId::Detector thedet = scRef->seed()->hitsAndFractions()[0].first.det();
1066 int subdet = scRef->seed()->hitsAndFractions()[0].first.subdetId();
1067 if (subdet == EcalBarrel) {
1068 preselCutValues = preselCutValuesBarrel_;
1069 } else if (subdet == EcalEndcap) {
1070 preselCutValues = preselCutValuesEndcap_;
1071 } else if (EcalTools::isHGCalDet(thedet)) {
1072 preselCutValues = preselCutValuesEndcap_;
1073 } else {
1074 edm::LogWarning("") << "GEDPhotonProducer: do not know if it is a barrel or endcap SuperCluster" << thedet << ' '
1075 << subdet;
1076 }
1077
1078
1079 if (parentSCRef.isNonnull() &&
1080 ptFast(parentSCRef->energy(), parentSCRef->position(), {0, 0, 0}) <= preselCutValues[0])
1081 continue;
1082
1083 reco::Photon newCandidate(*phoRef);
1084 iSC++;
1085
1086 if (runMVABasedHaloTagger_) {
1087 float BHmva = photonMVABasedHaloTagger_->calculateMVA(&newCandidate, globalCache()->haloTaggerGBR.get(), evt, es);
1088 newCandidate.setHaloTaggerMVAVal(BHmva);
1089 }
1090
1091
1092 reco::Photon::PflowIsolationVariables pfIso;
1093
1094
1095
1096 edm::Ptr<reco::Photon> photonPtr(photonHandle, lSC);
1097
1098 if (!recoStep_.isOOT()) {
1099 pfIso.chargedHadronIso = (*chargedHadrons)[photonPtr];
1100 pfIso.neutralHadronIso = (*neutralHadrons)[photonPtr];
1101 pfIso.photonIso = (*photons)[photonPtr];
1102 pfIso.chargedHadronWorstVtxIso = (*chargedHadronsWorstVtx)[photonPtr];
1103 pfIso.chargedHadronWorstVtxGeomVetoIso = (*chargedHadronsWorstVtxGeomVeto)[photonPtr];
1104 pfIso.chargedHadronPFPVIso = (*chargedHadronsPFPV)[photonPtr];
1105 }
1106
1107
1108
1109 pfIso.sumEcalClusterEt = !phoPFECALClusIsolationToken_.isUninitialized() ? (*pfEcalClusters)[photonPtr] : 0.;
1110 pfIso.sumHcalClusterEt = !phoPFHCALClusIsolationToken_.isUninitialized() ? (*pfHcalClusters)[photonPtr] : 0.;
1111 newCandidate.setPflowIsolationVariables(pfIso);
1112
1113
1114 photonEnergyCorrector_->calculate(evt, newCandidate, subdet, vertexCollection, es);
1115 if (candidateP4type_ == "fromEcalEnergy") {
1116 newCandidate.setP4(newCandidate.p4(reco::Photon::ecal_photons));
1117 newCandidate.setCandidateP4type(reco::Photon::ecal_photons);
1118 newCandidate.setMass(0.0);
1119 } else if (candidateP4type_ == "fromRegression1") {
1120 newCandidate.setP4(newCandidate.p4(reco::Photon::regression1));
1121 newCandidate.setCandidateP4type(reco::Photon::regression1);
1122 newCandidate.setMass(0.0);
1123 } else if (candidateP4type_ == "fromRegression2") {
1124 newCandidate.setP4(newCandidate.p4(reco::Photon::regression2));
1125 newCandidate.setCandidateP4type(reco::Photon::regression2);
1126 newCandidate.setMass(0.0);
1127 } else if (candidateP4type_ == "fromRefinedSCRegression") {
1128 newCandidate.setP4(newCandidate.p4(reco::Photon::regression2));
1129 newCandidate.setCandidateP4type(reco::Photon::regression2);
1130 newCandidate.setMass(0.0);
1131 }
1132
1133 outputPhotonCollection.push_back(newCandidate);
1134 }
1135 }