Back to home page

Project CMSSW displayed by LXR
 
 

    

File indexing completed on 2024-11-25 02:29:55 

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