Back to home page

Project CMSSW displayed by LXR
 
 

    

File indexing completed on 2022-06-16 01:50:15 

0001 /** \class GEDPhotonProducer
0002  **  
0003  **
0004  **  \author Nancy Marinelli, U. of Notre Dame, US
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     // Here we will have to load the DNN PFID if present in the config
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     ///for MVA based beam halo tagger in the EE
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   // std::string PhotonCoreCollection_;
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   //for isolation with map-based veto
0154   edm::EDGetTokenT<edm::ValueMap<std::vector<reco::PFCandidateRef>>> particleBasedIsolationToken;
0155   //photon isolation sums
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   //AA
0173   //Flags and severities to be excluded from calculations
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   //MIP
0197   std::unique_ptr<PhotonMIPHaloTagger> photonMIPHaloTagger_ = nullptr;
0198   //MVA based Halo tagger for the EE photons
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   // additional configuration and helpers
0212   std::unique_ptr<ElectronHcalHelper> hcalHelperCone_;
0213   std::unique_ptr<ElectronHcalHelper> hcalHelperBc_;
0214   bool hcalRun2EffDepth_;
0215 
0216   // DNN for PFID photon enabled
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 }  // namespace
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     //OOT photons in legacy 80X re-miniAOD do not have PF cluster embeded into the reco object
0299     //to preserve 80X behaviour
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   // R9 value to decide converted/unconverted
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   //AA
0341   //Flags and Severities to be excluded from photon calculations
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   // cut values for pre-selection
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   //moved from beginRun to here, I dont see how this could cause harm as its just reading in the exactly same parameters each run
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   ///Get the set for PF cluster isolation calculator
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   // Register the product
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   // this method is supposed to create, initialize and return a CacheData instance
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   // Get the PhotonCore collection
0492   bool validPhotonCoreHandle = false;
0493   Handle<reco::PhotonCoreCollection> photonCoreHandle;
0494   bool validPhotonHandle = false;
0495   Handle<reco::PhotonCollection> photonHandle;
0496   //value maps for isolation
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     //get isolation objects
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     //OOT photons in legacy 80X re-miniAOD workflow dont have cluster isolation embed in them
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   // Get EcalRecHits
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   // Get the  PF refined cluster  collection
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     // Get the  PF candidates collection
0558     theEvent.getByToken(pfCandidates_, pfCandidateHandle);
0559     //OOT photons have no PF candidates so its not an error in this case
0560     if (!pfCandidateHandle.isValid() && !recoStep_.isOOT()) {
0561       throw cms::Exception("GEDPhotonProducer") << "Error! Can't get the pfCandidates";
0562     }
0563   }
0564 
0565   // get the geometry from the event setup:
0566   caloGeom_ = &eventSetup.getData(caloGeometryToken_);
0567 
0568   // prepare access to hcal data
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   // Get the primary event vertex
0578   const reco::VertexCollection dummyVC;
0579   auto const& vertexCollection{usePrimaryVertex_ ? theEvent.get(vertexProducer_) : dummyVC};
0580 
0581   //  math::XYZPoint vtx(0.,0.,0.);
0582   //if (vertexCollection.size()>0) vtx = vertexCollection.begin()->position();
0583 
0584   // get the regression calculator ready
0585   photonEnergyCorrector_->init(eventSetup);
0586   if (photonEnergyCorrector_->gedRegression()) {
0587     photonEnergyCorrector_->gedRegression()->setEvent(theEvent);
0588     photonEnergyCorrector_->gedRegression()->setEventContent(eventSetup);
0589   }
0590 
0591   ///PF ECAL cluster based isolations
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;  // index in photon collection
0610   // Loop over barrel and endcap SC collections and fill the  photon collection
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                          //vtx,
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   // put the product in the event
0647   edm::LogInfo("GEDPhotonProducer") << " Put in the event " << iSC << " Photon Candidates \n";
0648 
0649   // go back to run2-like 2 effective depths if desired - depth 1 is the normal depth 1, depth 2 is the sum over the rest
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     //// Define the value map which associate to each  Egamma-unbiassaed candidate (key-ref) the corresponding PhotonRef
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     //// Fill the value map which associate to each Photon (key) the corresponding Egamma-unbiassaed candidate (value-ref)
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     //    const reco::SuperCluster* pClus=&(*scRef);
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     // SC energy preselection
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     //AA
0737     //Change these to consider severity level of hits
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     // fractional local covariances
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     //AA
0753     //Change these to consider severity level of hits
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     // for full5x5 local covariances, do noise-cleaning
0762     // by passing per crystal PF recHit thresholds and mult values.
0763     // mult values for EB and EE were obtained by dedicated studies.
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     // compute position of ECAL shower
0779     math::XYZPoint caloPosition = scRef->position();
0780 
0781     //// energy determination -- Default to create the candidate. Afterwards corrections are applied
0782     double photonEnergy = 1.;
0783     math::XYZPoint vtx(0., 0., 0.);
0784     if (!vertexCollection.empty())
0785       vtx = vertexCollection.begin()->position();
0786     // compute momentum vector of photon from primary vertex and cluster position
0787     math::XYZVector direction = caloPosition - vtx;
0788     //math::XYZVector momentum = direction.unit() * photonEnergy ;
0789     math::XYZVector momentum = direction.unit();
0790 
0791     // Create dummy candidate with unit momentum and zero energy to allow setting of all variables. The energy is set for last.
0792     math::XYZTLorentzVectorD p4(momentum.x(), momentum.y(), momentum.z(), photonEnergy);
0793     reco::Photon newCandidate(p4, caloPosition, coreRef, vtx);
0794 
0795     //std::cout << " standard p4 before " << newCandidate.p4() << " energy " << newCandidate.energy() <<  std::endl;
0796     //std::cout << " type " <<newCandidate.getCandidateP4type() <<  " standard p4 after " << newCandidate.p4() << " energy " << newCandidate.energy() << std::endl;
0797 
0798     // Calculate fiducial flags and isolation variable. Blocked are filled from the isolationCalculator
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     /// fill shower shape block
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     /// fill extra shower shapes
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     // fill preshower shapes
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     /// fill full5x5 shower shape block
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     /// fill extra full5x5 shower shapes
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     // fill preshower shapes
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     //get the pointer for the photon object
0941     edm::Ptr<reco::PhotonCore> photonPtr(photonCoreHandle, lSC);
0942 
0943     // New in CMSSW_12_1_0 for PFID with DNNs
0944     // The PFIso values are computed in the first loop on gedPhotonsTmp to make them available as DNN inputs.
0945     // They are computed with the same inputs and algo as the final PFiso variables computed in the second loop after PF.
0946     // Get PFClusters for PFID only if the PFID DNN evaluation is enabled
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     /// get ecal photon specific corrected energy
0962     /// plus values from regressions     and store them in the Photon
0963     // Photon candidate takes by default (set in photons_cfi.py)
0964     // a 4-momentum derived from the ecal photon-specific corrections.
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       // Make it an EE photon
0990       reco::Photon::FiducialFlags fiducialFlags;
0991       fiducialFlags.isEE = true;
0992       newCandidate.setFiducialVolumeFlags(fiducialFlags);
0993     }
0994 
0995     // fill MIP Vairables for Halo: Block for MIP are filled from PhotonMIPHaloTagger
0996     reco::Photon::MIPVariables mipVar;
0997     if (subdet == EcalBarrel && runMIPTagger_) {
0998       photonMIPHaloTagger_->MIPcalculate(&newCandidate, evt, es, mipVar);
0999       newCandidate.setMIPVariables(mipVar);
1000     }
1001 
1002     /// Pre-selection loose  isolation cuts
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     // Here send the list of photons to the PhotonDNNEstimator and get back the values for all the photons in one go
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       // The model index it is not useful for the moment
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     // SC energy preselection
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_) {  ///sets values only for EE, for EB it always returns 1
1087       float BHmva = photonMVABasedHaloTagger_->calculateMVA(&newCandidate, globalCache()->haloTaggerGBR.get(), evt, es);
1088       newCandidate.setHaloTaggerMVAVal(BHmva);
1089     }
1090 
1091     // Calculate the PF isolation
1092     reco::Photon::PflowIsolationVariables pfIso;
1093     // The PFID are not recomputed since they have been already computed in the first loop with the DNN
1094 
1095     //get the pointer for the photon object
1096     edm::Ptr<reco::Photon> photonPtr(photonHandle, lSC);
1097 
1098     if (!recoStep_.isOOT()) {  //out of time photons do not have PF info so skip in this case
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     //OOT photons in legacy 80X reminiAOD workflow dont have pf cluster isolation embeded into them at this stage
1108     // They have been already computed in the first loop on gedPhotonsTmp but better to compute them again here.
1109     pfIso.sumEcalClusterEt = !phoPFECALClusIsolationToken_.isUninitialized() ? (*pfEcalClusters)[photonPtr] : 0.;
1110     pfIso.sumHcalClusterEt = !phoPFHCALClusIsolationToken_.isUninitialized() ? (*pfHcalClusters)[photonPtr] : 0.;
1111     newCandidate.setPflowIsolationVariables(pfIso);
1112 
1113     // do the regression
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 }
[ Source navigation ] [ Diff markup ] [ Identifier search ] [ general search ]

This page was automatically generated by the 2.2.1 LXR engine.

The LXR team
Valid CSS 2.1! Valid HTML 4.01!