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