Back to home page

Project CMSSW displayed by LXR
 
 

    

File indexing completed on 2023-11-29 01:33:07 

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