Implementation of beam halo tagger algorithm for the endcap photons

DataFormats/EgammaCandidates/interface/Photon.h

@@ -579,6 +579,12 @@ namespace reco {
     // 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
     void hcalToRun2EffDepth();
 
+    ///MVA based beam halo tagger - trained for EE and for pT > 200 GeV
+    float haloTaggerMVAVal() const { return haloTaggerMVAVal_; }
+
+    ///set the haloTaggerMVAVal here
+    void setHaloTaggerMVAVal(float x) { haloTaggerMVAVal_ = x; }
+
   private:
     /// check overlap with another candidate
     bool overlap(const Candidate&) const override;
@@ -599,6 +605,7 @@ namespace reco {
     MIPVariables mipVariableBlock_;
     PflowIsolationVariables pfIsolation_;
     PflowIDVariables pfID_;
+    float haloTaggerMVAVal_;
   };
 
 }  // namespace reco

DataFormats/EgammaCandidates/src/Photon.cc

@@ -4,7 +4,11 @@
 using namespace reco;
 
 Photon::Photon(const LorentzVector& p4, const Point& caloPos, const PhotonCoreRef& core, const Point& vtx)
-    : RecoCandidate(0, p4, vtx, 22), caloPosition_(caloPos), photonCore_(core), pixelSeed_(false) {}
+    : RecoCandidate(0, p4, vtx, 22),
+      caloPosition_(caloPos),
+      photonCore_(core),
+      pixelSeed_(false),
+      haloTaggerMVAVal_(-999) {}
 
 Photon::Photon(const Photon& rhs)
     : RecoCandidate(rhs),
@@ -19,7 +23,8 @@ Photon::Photon(const Photon& rhs)
       saturationInfo_(rhs.saturationInfo_),
       eCorrections_(rhs.eCorrections_),
       mipVariableBlock_(rhs.mipVariableBlock_),
-      pfIsolation_(rhs.pfIsolation_) {}
+      pfIsolation_(rhs.pfIsolation_),
+      haloTaggerMVAVal_(rhs.haloTaggerMVAVal_) {}
 
 Photon::~Photon() {}
 

DataFormats/EgammaCandidates/src/classes_def.xml

@@ -9,7 +9,8 @@
   </ioread>
 
 
-  <class name="reco::Photon" ClassVersion="17">
+  <class name="reco::Photon" ClassVersion="18">
+   <version ClassVersion="18" checksum="2393532952"/>
    <version ClassVersion="17" checksum="4077732720"/>
    <version ClassVersion="16" checksum="3459546220"/>
    <version ClassVersion="15" checksum="577991108"/>

DataFormats/PatCandidates/src/classes_def_objects.xml

@@ -222,7 +222,8 @@
   </ioread>
 
   <class name="std::vector<pat::tau::TauPFEssential>" />
-  <class name="pat::Photon"  ClassVersion="21">
+  <class name="pat::Photon"  ClassVersion="22">
+   <version ClassVersion="22" checksum="1415449268"/>
    <version ClassVersion="21" checksum="3263223164"/>
    <version ClassVersion="20" checksum="1579185367"/>
    <version ClassVersion="19" checksum="4285818507"/>

RecoEgamma/EgammaPhotonProducers/python/gedPhotons_cfi.py

@@ -3,6 +3,7 @@
 from RecoEgamma.PhotonIdentification.pfIsolationCalculator_cfi import *
 from RecoEgamma.PhotonIdentification.isolationCalculator_cfi import *
 from RecoEgamma.PhotonIdentification.mipVariable_cfi import *
+from RecoEgamma.PhotonIdentification.mvaHaloVariable_cfi import *
 from RecoEcal.EgammaClusterProducers.hybridSuperClusters_cfi import *
 from RecoEcal.EgammaClusterProducers.multi5x5BasicClusters_cfi import *
 
@@ -39,6 +40,7 @@
     isolationSumsCalculatorSet = cms.PSet(isolationSumsCalculator),
     PFIsolationCalculatorSet = cms.PSet(pfIsolationCalculator),                        
     mipVariableSet = cms.PSet(mipVariable), 
+    mvaBasedHaloVariableSet = cms.PSet(mvaHaloVariable), 
     usePrimaryVertex = cms.bool(True),
     primaryVertexProducer = cms.InputTag('offlinePrimaryVerticesWithBS'),
     posCalc_t0_endcPresh = cms.double(3.6),
@@ -51,6 +53,7 @@
     endcapEcalHits = cms.InputTag("ecalRecHit","EcalRecHitsEE"),
     preshowerHits = cms.InputTag("ecalPreshowerRecHit","EcalRecHitsES"),
     runMIPTagger = cms.bool(True),
+    runMVABasedHaloTagger = cms.bool(True),
     highEt  = cms.double(100.),                       
     minR9Barrel = cms.double(0.94),
     minR9Endcap = cms.double(0.95),

RecoEgamma/EgammaPhotonProducers/python/photons_cfi.py

@@ -2,6 +2,7 @@
 
 from RecoEgamma.PhotonIdentification.isolationCalculator_cfi import *
 from RecoEgamma.PhotonIdentification.mipVariable_cfi import *
+from RecoEgamma.PhotonIdentification.mvaHaloVariable_cfi import *
 from RecoEcal.EgammaClusterProducers.hybridSuperClusters_cfi import *
 from RecoEcal.EgammaClusterProducers.multi5x5BasicClusters_cfi import *
 from RecoEgamma.EgammaIsolationAlgos.egammaHBHERecHitThreshold_cff import egammaHBHERecHit
@@ -30,6 +31,7 @@
     candidateP4type = cms.string("fromEcalEnergy"),                     
     isolationSumsCalculatorSet = cms.PSet(isolationSumsCalculator),
     mipVariableSet = cms.PSet(mipVariable), 
+    mvaBasedHaloVariableSet = cms.PSet(mvaHaloVariable), 
     usePrimaryVertex = cms.bool(True),
     primaryVertexProducer = cms.InputTag('offlinePrimaryVerticesWithBS'),
     posCalc_t0_endcPresh = cms.double(3.6),
@@ -42,6 +44,7 @@
     endcapEcalHits = cms.InputTag("ecalRecHit","EcalRecHitsEE"),
     preshowerHits = cms.InputTag("ecalPreshowerRecHit","EcalRecHitsES"),
     runMIPTagger = cms.bool(True),
+    runMVABasedHaloTagger = cms.bool(False),
     highEt  = cms.double(100.),                       
     minR9Barrel = cms.double(0.94),
     minR9Endcap = cms.double(0.95),
@@ -161,6 +164,7 @@
     hbheModule = cms.string('hbhereco'),
     endcapEcalHits = cms.InputTag("ecalRecHit","EcalRecHitsEE"),
     runMIPTagger = cms.bool(True),
+    runMVABasedHaloTagger = cms.bool(False),
     highEt  = cms.double(100.),
     minR9Barrel = cms.double(10.0),
     minR9Endcap = cms.double(10.0),

RecoEgamma/EgammaPhotonProducers/src/GEDPhotonProducer.cc

@@ -40,6 +40,7 @@
 #include "RecoEgamma/EgammaPhotonAlgos/interface/PhotonEnergyCorrector.h"
 #include "RecoEgamma/PhotonIdentification/interface/PhotonIsolationCalculator.h"
 #include "RecoEgamma/PhotonIdentification/interface/PhotonMIPHaloTagger.h"
+#include "RecoEgamma/PhotonIdentification/interface/PhotonMVABasedHaloTagger.h"
 #include "RecoLocalCalo/EcalRecAlgos/interface/EcalSeverityLevelAlgo.h"
 #include "RecoLocalCalo/EcalRecAlgos/interface/EcalSeverityLevelAlgoRcd.h"
 #include "CondFormats/EcalObjects/interface/EcalPFRecHitThresholds.h"
@@ -50,6 +51,8 @@
 #include "PhysicsTools/TensorFlow/interface/TensorFlow.h"
 #include "RecoEgamma/EgammaIsolationAlgos/interface/EcalPFClusterIsolation.h"
 #include "RecoEgamma/EgammaIsolationAlgos/interface/HcalPFClusterIsolation.h"
+#include "CondFormats/GBRForest/interface/GBRForest.h"
+#include "CommonTools/MVAUtils/interface/GBRForestTools.h"
 
 class CacheData {
 public:
@@ -67,8 +70,16 @@ class CacheData {
       const auto useEBModelInGap = pset_dnn.getParameter<bool>("useEBModelInGap");
       photonDNNEstimator = std::make_unique<PhotonDNNEstimator>(config, useEBModelInGap);
     }
+    ///for MVA based beam halo tagger in the EE
+    const auto runMVABasedHaloTagger = conf.getParameter<bool>("runMVABasedHaloTagger");
+    edm::ParameterSet mvaBasedHaloVariableSet = conf.getParameter<edm::ParameterSet>("mvaBasedHaloVariableSet");
+    auto trainingFileName_ = mvaBasedHaloVariableSet.getParameter<edm::FileInPath>(("trainingFileName")).fullPath();
+    if (runMVABasedHaloTagger) {
+      haloTaggerGBR = createGBRForest(trainingFileName_);
+    }
   }
   std::unique_ptr<const PhotonDNNEstimator> photonDNNEstimator;
+  std::unique_ptr<const GBRForest> haloTaggerGBR;
 };
 
 class GEDPhotonProducer : public edm::stream::EDProducer<edm::GlobalCache<CacheData>> {
@@ -174,6 +185,7 @@ class GEDPhotonProducer : public edm::stream::EDProducer<edm::GlobalCache<CacheD
   double minR9Barrel_;
   double minR9Endcap_;
   bool runMIPTagger_;
+  bool runMVABasedHaloTagger_;
 
   RecoStepInfo recoStep_;
 
@@ -183,6 +195,8 @@ class GEDPhotonProducer : public edm::stream::EDProducer<edm::GlobalCache<CacheD
 
   //MIP
   std::unique_ptr<PhotonMIPHaloTagger> photonMIPHaloTagger_ = nullptr;
+  //MVA based Halo tagger for the EE photons
+  std::unique_ptr<PhotonMVABasedHaloTagger> photonMVABasedHaloTagger_ = nullptr;
 
   std::vector<double> preselCutValuesBarrel_;
   std::vector<double> preselCutValuesEndcap_;
@@ -318,6 +332,7 @@ GEDPhotonProducer::GEDPhotonProducer(const edm::ParameterSet& config, const Cach
   minR9Endcap_ = config.getParameter<double>("minR9Endcap");
   usePrimaryVertex_ = config.getParameter<bool>("usePrimaryVertex");
   runMIPTagger_ = config.getParameter<bool>("runMIPTagger");
+  runMVABasedHaloTagger_ = config.getParameter<bool>("runMVABasedHaloTagger");
 
   candidateP4type_ = config.getParameter<std::string>("candidateP4type");
   valueMapPFCandPhoton_ = config.getParameter<std::string>("valueMapPhotons");
@@ -411,6 +426,12 @@ GEDPhotonProducer::GEDPhotonProducer(const edm::ParameterSet& config, const Cach
     photonMIPHaloTagger_->setup(mipVariableSet, consumesCollector());
   }
 
+  if (recoStep_.isFinal() && runMVABasedHaloTagger_) {
+    edm::ParameterSet mvaBasedHaloVariableSet = config.getParameter<edm::ParameterSet>("mvaBasedHaloVariableSet");
+    photonMVABasedHaloTagger_ =
+        std::make_unique<PhotonMVABasedHaloTagger>(mvaBasedHaloVariableSet, consumesCollector());
+  }
+
   ///Get the set for PF cluster isolation calculator
   const edm::ParameterSet& pfECALClusIsolCfg = config.getParameter<edm::ParameterSet>("pfECALClusIsolCfg");
   pfClusterProducer_ =
@@ -1057,6 +1078,11 @@ void GEDPhotonProducer::fillPhotonCollection(edm::Event& evt,
     reco::Photon newCandidate(*phoRef);
     iSC++;
 
+    if (runMVABasedHaloTagger_) {  ///sets values only for EE, for EB it always returns 1
+      float BHmva = photonMVABasedHaloTagger_->calculateMVA(&newCandidate, globalCache()->haloTaggerGBR.get(), evt, es);
+      newCandidate.setHaloTaggerMVAVal(BHmva);
+    }
+
     // Calculate the PF isolation
     reco::Photon::PflowIsolationVariables pfIso;
     // The PFID are not recomputed since they have been already computed in the first loop with the DNN

RecoEgamma/EgammaPhotonProducers/src/PhotonProducer.cc

@@ -34,6 +34,7 @@
 #include "RecoEgamma/EgammaPhotonAlgos/interface/PhotonEnergyCorrector.h"
 #include "RecoEgamma/PhotonIdentification/interface/PhotonIsolationCalculator.h"
 #include "RecoEgamma/PhotonIdentification/interface/PhotonMIPHaloTagger.h"
+#include "RecoEgamma/PhotonIdentification/interface/PhotonMVABasedHaloTagger.h"
 #include "RecoLocalCalo/EcalRecAlgos/interface/EcalSeverityLevelAlgo.h"
 #include "RecoLocalCalo/EcalRecAlgos/interface/EcalSeverityLevelAlgoRcd.h"
 #include "RecoEgamma/EgammaElectronAlgos/interface/ElectronHcalHelper.h"
@@ -83,6 +84,7 @@ class PhotonProducer : public edm::stream::EDProducer<> {
   double minR9Barrel_;
   double minR9Endcap_;
   bool runMIPTagger_;
+  bool runMVABasedHaloTagger_;
 
   bool validConversions_;
 
@@ -95,6 +97,8 @@ class PhotonProducer : public edm::stream::EDProducer<> {
 
   //MIP
   PhotonMIPHaloTagger photonMIPHaloTagger_;
+  //MVA based Halo tagger for the EE photons
+  std::unique_ptr<PhotonMVABasedHaloTagger> photonMVABasedHaloTagger_ = nullptr;
 
   std::vector<double> preselCutValuesBarrel_;
   std::vector<double> preselCutValuesEndcap_;
@@ -127,6 +131,7 @@ PhotonProducer::PhotonProducer(const edm::ParameterSet& config)
   minR9Endcap_ = config.getParameter<double>("minR9Endcap");
   usePrimaryVertex_ = config.getParameter<bool>("usePrimaryVertex");
   runMIPTagger_ = config.getParameter<bool>("runMIPTagger");
+  runMVABasedHaloTagger_ = config.getParameter<bool>("runMVABasedHaloTagger");
 
   candidateP4type_ = config.getParameter<std::string>("candidateP4type");
 

RecoEgamma/PhotonIdentification/BuildFile.xml

@@ -24,6 +24,7 @@
 <use name="RecoLocalCalo/EcalRecAlgos"/>
 <use name="RecoLocalCalo/HcalRecAlgos"/>
 <use name="PhysicsTools/TensorFlow" />
+<use name="CommonTools/MVAUtils" />
 <export>
   <lib name="1"/>
 </export>

RecoEgamma/PhotonIdentification/interface/PhotonMVABasedHaloTagger.h

@@ -0,0 +1,90 @@
+/** \class PhotonMVABasedHaloTagger
+ *   \author Shilpi Jain (University of Minnesota)
+ * Links to the presentation: 
+ 1. ECAL DPG: https://indico.cern.ch/event/991261/contributions/4283096/attachments/2219229/3757719/beamHalo_31march_v1.pdf
+ 2. JetMET POG: https://indico.cern.ch/event/1027614/contributions/4314949/attachments/2224472/3767396/beamHalo_12April.pdf
+ */
+
+#ifndef RecoEgamma_PhotonIdentification_PhotonMVABasedHaloTagger_h
+#define RecoEgamma_PhotonIdentification_PhotonMVABasedHaloTagger_h
+
+#include "FWCore/Framework/interface/Event.h"
+#include "FWCore/Framework/interface/EventSetup.h"
+#include "FWCore/ParameterSet/interface/ParameterSet.h"
+#include "DataFormats/EcalRecHit/interface/EcalRecHitCollections.h"
+#include "DataFormats/EgammaCandidates/interface/Photon.h"
+#include "FWCore/Framework/interface/ConsumesCollector.h"
+#include "Geometry/CaloGeometry/interface/CaloGeometry.h"
+#include "Geometry/Records/interface/CaloGeometryRecord.h"
+#include "DataFormats/HcalRecHit/interface/HcalRecHitCollections.h"
+#include "CondFormats/GBRForest/interface/GBRForest.h"
+#include "RecoEcal/EgammaCoreTools/interface/EcalClusterLazyTools.h"
+#include "RecoEgamma/EgammaIsolationAlgos/interface/EgammaHcalIsolation.h"
+#include "CondFormats/EcalObjects/interface/EcalPFRecHitThresholds.h"
+#include "CondFormats/DataRecord/interface/EcalPFRecHitThresholdsRcd.h"
+#include <vector>
+
+class PhotonMVABasedHaloTagger {
+public:
+  PhotonMVABasedHaloTagger(const edm::ParameterSet& conf, edm::ConsumesCollector&& iC);
+
+  double calculateMVA(const reco::Photon* pho,
+                      const GBRForest* gbr_,
+                      const edm::Event& iEvent,
+                      const edm::EventSetup& es);
+
+private:
+  void calphoClusCoordinECAL(const CaloGeometry* geo,
+                             const reco::Photon*,
+                             const EcalPFRecHitThresholds* thresholds,
+                             const EcalRecHitCollection& ecalRecHits);
+
+  void calmatchedHBHECoordForBothHypothesis(const CaloGeometry* geo,
+                                            const reco::Photon*,
+                                            const HBHERecHitCollection& HBHERecHits);
+
+  void calmatchedESCoordForBothHypothesis(const CaloGeometry* geo,
+                                          const reco::Photon*,
+                                          const EcalRecHitCollection& ESRecHits);
+
+  double calAngleBetweenEEAndSubDet(int nhits, double subdetClusX, double subdetClusY, double subdetClusZ);
+
+  int hcalClusNhits_samedPhi_, hcalClusNhits_samedR_;
+  int ecalClusNhits_, preshowerNhits_samedPhi_, preshowerNhits_samedR_;
+  double hcalClusX_samedPhi_, hcalClusY_samedPhi_, hcalClusZ_samedPhi_, hcalClusX_samedR_, hcalClusY_samedR_,
+      hcalClusZ_samedR_;
+  double hcalClusE_samedPhi_, hcalClusE_samedR_;
+
+  double ecalClusX_, ecalClusY_, ecalClusZ_;
+  double preshowerX_samedPhi_, preshowerY_samedPhi_, preshowerZ_samedPhi_, preshowerX_samedR_, preshowerY_samedR_,
+      preshowerZ_samedR_;
+  double ecalClusE_, preshowerE_samedPhi_, preshowerE_samedR_;
+  double noiseThrES_;
+
+  EgammaHcalIsolation::arrayHB recHitEThresholdHB_;
+  EgammaHcalIsolation::arrayHE recHitEThresholdHE_;
+
+  const edm::ESGetToken<CaloGeometry, CaloGeometryRecord> geometryToken_;
+  const edm::ESGetToken<EcalPFRecHitThresholds, EcalPFRecHitThresholdsRcd> ecalPFRechitThresholdsToken_;
+  const EcalClusterLazyTools::ESGetTokens ecalClusterToolsESGetTokens_;
+
+  edm::ESHandle<CaloGeometry> pG_;
+  edm::EDGetTokenT<double> rhoLabel_;
+  edm::EDGetTokenT<EcalRecHitCollection> EBecalCollection_;
+  edm::EDGetTokenT<EcalRecHitCollection> EEecalCollection_;
+  edm::EDGetTokenT<EcalRecHitCollection> ESCollection_;
+  edm::EDGetTokenT<HBHERecHitCollection> HBHERecHitsCollection_;
+
+  ///values of dR etc to cluster the hits in various sub-detectors
+  static constexpr float dr2Max_ECALClus_ = 0.2 * 0.2;
+  static constexpr float rho2Min_ECALpos_ = 31 * 31;            //cm
+  static constexpr float rho2Max_ECALpos_ = 172 * 172;          //cm
+  static constexpr float dRho2Max_HCALClus_SamePhi_ = 26 * 26;  //cm
+  static constexpr float dPhiMax_HCALClus_SamePhi_ = 0.15;
+  static constexpr float dR2Max_HCALClus_SamePhi_ = 0.15 * 0.15;
+  static constexpr float dRho2Max_ESClus_ = 2.2 * 2.2;  //cm
+  static constexpr float dXY_ESClus_SamePhi_ = 1;       ///cm
+  static constexpr float dXY_ESClus_SamedR_ = 1;        ///cm
+};
+
+#endif  // PhotonMVABasedHaloTagger_H

RecoEgamma/PhotonIdentification/python/mvaHaloVariable_cfi.py

@@ -0,0 +1,20 @@
+import FWCore.ParameterSet.Config as cms
+from RecoEgamma.EgammaIsolationAlgos.egammaHBHERecHitThreshold_cff import egammaHBHERecHit
+
+pathToHaloMVATrainingFile = "RecoEgamma/PhotonIdentification/data/beamHaloTaggerID/xgboostToTMVA_BHtagger.root"
+mvaHaloVariable = cms.PSet(
+    #required inputs
+    trainingFileName = cms.FileInPath(pathToHaloMVATrainingFile),
+    rhoLabel = cms.InputTag("fixedGridRhoFastjetAllTmp"),
+    barrelEcalRecHitCollection = cms.InputTag("ecalRecHit","EcalRecHitsEB"),
+    endcapEcalRecHitCollection = cms.InputTag("ecalRecHit","EcalRecHitsEE"),
+    esRecHitCollection = cms.InputTag("ecalPreshowerRecHit","EcalRecHitsES"),
+    HBHERecHitsCollection = egammaHBHERecHit.hbheRecHits,
+    recHitEThresholdHB = egammaHBHERecHit.recHitEThresholdHB,
+    recHitEThresholdHE = egammaHBHERecHit.recHitEThresholdHE,
+    noiseThrES = cms.double(0.0)
+
+
+)
+
+