From b143f4b0774bc5830c89a7bbbb9df7a9ad8b3724 Mon Sep 17 00:00:00 2001 From: Shilpi Date: Wed, 2 Feb 2022 11:01:58 +0100 Subject: [PATCH 01/13] MVA based EE photon tagger implementation --- .../EgammaCandidates/interface/Photon.h | 9 + DataFormats/EgammaCandidates/src/Photon.cc | 5 +- .../EgammaCandidates/src/classes_def.xml | 3 +- .../PatCandidates/src/classes_def_objects.xml | 3 +- .../python/gedPhotons_cfi.py | 3 + .../python/photons_cfi.py | 3 + .../src/GEDPhotonProducer.cc | 16 + .../src/PhotonProducer.cc | 16 + RecoEgamma/PhotonIdentification/BuildFile.xml | 1 + .../interface/PhotonMVABasedHaloTagger.h | 92 ++++ .../python/mvaHaloVariable_cfi.py | 20 + .../src/PhotonMVABasedHaloTagger.cc | 468 ++++++++++++++++++ 12 files changed, 635 insertions(+), 4 deletions(-) create mode 100644 RecoEgamma/PhotonIdentification/interface/PhotonMVABasedHaloTagger.h create mode 100644 RecoEgamma/PhotonIdentification/python/mvaHaloVariable_cfi.py create mode 100644 RecoEgamma/PhotonIdentification/src/PhotonMVABasedHaloTagger.cc diff --git a/DataFormats/EgammaCandidates/interface/Photon.h b/DataFormats/EgammaCandidates/interface/Photon.h index f1daf95615eb4..f4b3bc4e2bc1c 100644 --- a/DataFormats/EgammaCandidates/interface/Photon.h +++ b/DataFormats/EgammaCandidates/interface/Photon.h @@ -579,6 +579,14 @@ 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 +607,7 @@ namespace reco { MIPVariables mipVariableBlock_; PflowIsolationVariables pfIsolation_; PflowIDVariables pfID_; + float haloTaggerMVAVal_; }; } // namespace reco diff --git a/DataFormats/EgammaCandidates/src/Photon.cc b/DataFormats/EgammaCandidates/src/Photon.cc index f1ac9770a17fe..27892904b9144 100644 --- a/DataFormats/EgammaCandidates/src/Photon.cc +++ b/DataFormats/EgammaCandidates/src/Photon.cc @@ -4,7 +4,7 @@ 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 +19,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() {} diff --git a/DataFormats/EgammaCandidates/src/classes_def.xml b/DataFormats/EgammaCandidates/src/classes_def.xml index 9432329160786..20137d71f63d4 100644 --- a/DataFormats/EgammaCandidates/src/classes_def.xml +++ b/DataFormats/EgammaCandidates/src/classes_def.xml @@ -9,7 +9,8 @@ - + + diff --git a/DataFormats/PatCandidates/src/classes_def_objects.xml b/DataFormats/PatCandidates/src/classes_def_objects.xml index 2b98aef02aab2..b4a6802919c71 100644 --- a/DataFormats/PatCandidates/src/classes_def_objects.xml +++ b/DataFormats/PatCandidates/src/classes_def_objects.xml @@ -222,7 +222,8 @@ - + + diff --git a/RecoEgamma/EgammaPhotonProducers/python/gedPhotons_cfi.py b/RecoEgamma/EgammaPhotonProducers/python/gedPhotons_cfi.py index d8d6f7b86e32f..cf6dfa5e5ebf5 100644 --- a/RecoEgamma/EgammaPhotonProducers/python/gedPhotons_cfi.py +++ b/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), diff --git a/RecoEgamma/EgammaPhotonProducers/python/photons_cfi.py b/RecoEgamma/EgammaPhotonProducers/python/photons_cfi.py index 1ffac6650827f..bc05e54af71c8 100644 --- a/RecoEgamma/EgammaPhotonProducers/python/photons_cfi.py +++ b/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(True), highEt = cms.double(100.), minR9Barrel = cms.double(0.94), minR9Endcap = cms.double(0.95), diff --git a/RecoEgamma/EgammaPhotonProducers/src/GEDPhotonProducer.cc b/RecoEgamma/EgammaPhotonProducers/src/GEDPhotonProducer.cc index 337cea3032416..30713677ace12 100644 --- a/RecoEgamma/EgammaPhotonProducers/src/GEDPhotonProducer.cc +++ b/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" @@ -174,6 +175,7 @@ class GEDPhotonProducer : public edm::stream::EDProducer photonMIPHaloTagger_ = nullptr; + //MVA based Halo tagger for the EE photons + std::unique_ptr photonMVABasedHaloTagger_ = nullptr; std::vector preselCutValuesBarrel_; std::vector preselCutValuesEndcap_; @@ -318,6 +322,7 @@ GEDPhotonProducer::GEDPhotonProducer(const edm::ParameterSet& config, const Cach minR9Endcap_ = config.getParameter("minR9Endcap"); usePrimaryVertex_ = config.getParameter("usePrimaryVertex"); runMIPTagger_ = config.getParameter("runMIPTagger"); + runMVABasedHaloTagger_ = config.getParameter("runMVABasedHaloTagger"); candidateP4type_ = config.getParameter("candidateP4type"); valueMapPFCandPhoton_ = config.getParameter("valueMapPhotons"); @@ -411,6 +416,12 @@ GEDPhotonProducer::GEDPhotonProducer(const edm::ParameterSet& config, const Cach photonMIPHaloTagger_->setup(mipVariableSet, consumesCollector()); } + if (recoStep_.isFinal() && runMVABasedHaloTagger_) { + edm::ParameterSet mvaBasedHaloVariableSet = config.getParameter("mvaBasedHaloVariableSet"); + photonMVABasedHaloTagger_ = std::make_unique(mvaBasedHaloVariableSet, consumesCollector()); + + } + ///Get the set for PF cluster isolation calculator const edm::ParameterSet& pfECALClusIsolCfg = config.getParameter("pfECALClusIsolCfg"); pfClusterProducer_ = @@ -1057,6 +1068,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, 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 diff --git a/RecoEgamma/EgammaPhotonProducers/src/PhotonProducer.cc b/RecoEgamma/EgammaPhotonProducers/src/PhotonProducer.cc index 8fdf6f536cbe2..4db6d3e9b0074 100644 --- a/RecoEgamma/EgammaPhotonProducers/src/PhotonProducer.cc +++ b/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_ = nullptr; std::vector preselCutValuesBarrel_; std::vector preselCutValuesEndcap_; @@ -127,6 +131,7 @@ PhotonProducer::PhotonProducer(const edm::ParameterSet& config) minR9Endcap_ = config.getParameter("minR9Endcap"); usePrimaryVertex_ = config.getParameter("usePrimaryVertex"); runMIPTagger_ = config.getParameter("runMIPTagger"); + runMVABasedHaloTagger_ = config.getParameter("runMVABasedHaloTagger"); candidateP4type_ = config.getParameter("candidateP4type"); @@ -231,6 +236,12 @@ PhotonProducer::PhotonProducer(const edm::ParameterSet& config) edm::ParameterSet mipVariableSet = config.getParameter("mipVariableSet"); photonMIPHaloTagger_.setup(mipVariableSet, consumesCollector()); + if (runMVABasedHaloTagger_) { + edm::ParameterSet mvaBasedHaloVariableSet = config.getParameter("mvaBasedHaloVariableSet"); + photonMVABasedHaloTagger_ = std::make_unique(mvaBasedHaloVariableSet, consumesCollector()); + } + + // Register the product produces(PhotonCollection_); } @@ -487,6 +498,11 @@ void PhotonProducer::fillPhotonCollection(edm::Event& evt, newCandidate.setMIPVariables(mipVar); } + if (runMVABasedHaloTagger_) { ///sets values only for EE, for EB it always returns 1 + double BHmva = photonMVABasedHaloTagger_->calculateMVA(&newCandidate, evt, es); + newCandidate.setHaloTaggerMVAVal(BHmva); + } + /// Pre-selection loose isolation cuts bool isLooseEM = true; if (newCandidate.pt() < highEt_) { diff --git a/RecoEgamma/PhotonIdentification/BuildFile.xml b/RecoEgamma/PhotonIdentification/BuildFile.xml index 02e53c5fd2adf..1aeb92479f3cf 100644 --- a/RecoEgamma/PhotonIdentification/BuildFile.xml +++ b/RecoEgamma/PhotonIdentification/BuildFile.xml @@ -24,6 +24,7 @@ + diff --git a/RecoEgamma/PhotonIdentification/interface/PhotonMVABasedHaloTagger.h b/RecoEgamma/PhotonIdentification/interface/PhotonMVABasedHaloTagger.h new file mode 100644 index 0000000000000..f3498f526fc8b --- /dev/null +++ b/RecoEgamma/PhotonIdentification/interface/PhotonMVABasedHaloTagger.h @@ -0,0 +1,92 @@ +/** \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 PhotonMVABasedHaloTagger_H +#define 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 + +class PhotonMVABasedHaloTagger { +public: + PhotonMVABasedHaloTagger(const edm::ParameterSet& conf, edm::ConsumesCollector&& iC); + virtual ~PhotonMVABasedHaloTagger(){}; + + double calculateMVA(const reco::Photon* pho, + const edm::Event&iEvent, + const edm::EventSetup& es + ); + + + void calphoClusCoordinECAL(const CaloGeometry* geo, + const reco::Photon*, + const EcalPFRecHitThresholds *thresholds, + edm::Handle ecalrechitCollHandle + ); + + + void calmatchedHBHECoordForBothHypothesis(const CaloGeometry* geo, + const reco::Photon*, + edm::Handle hbheHitHandle + ); + + + void calmatchedESCoordForBothHypothesis(const CaloGeometry* geo, + const reco::Photon*, + edm::Handle esrechitCollHandle + ); + + double calAngleBetweenEEAndSubDet(int nhits, + double subdetClusX, + double subdetClusY, + double subdetClusZ); + + +private: + 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_; + + std::string trainingFileName_; + + const edm::ESGetToken geometryToken_; + const edm::ESGetToken ecalPFRechitThresholdsToken_; + const EcalClusterLazyTools::ESGetTokens ecalClusterToolsESGetTokens_; + + edm::ESHandle pG_; + edm::EDGetTokenT rhoLabel_; + edm::EDGetTokenT EBecalCollection_; + edm::EDGetTokenT EEecalCollection_; + edm::EDGetTokenT ESCollection_; + edm::EDGetTokenT HBHERecHitsCollection_; + std::unique_ptr gbr_; + +}; + +#endif // PhotonMVABasedHaloTagger_H diff --git a/RecoEgamma/PhotonIdentification/python/mvaHaloVariable_cfi.py b/RecoEgamma/PhotonIdentification/python/mvaHaloVariable_cfi.py new file mode 100644 index 0000000000000..49b99cf575a9a --- /dev/null +++ b/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/xgboostToTMVA_BHtagger.xml" +mvaHaloVariable = cms.PSet( + #required inputs + trainingFileName = cms.string(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) + + +) + + diff --git a/RecoEgamma/PhotonIdentification/src/PhotonMVABasedHaloTagger.cc b/RecoEgamma/PhotonIdentification/src/PhotonMVABasedHaloTagger.cc new file mode 100644 index 0000000000000..9154c8a7efc9c --- /dev/null +++ b/RecoEgamma/PhotonIdentification/src/PhotonMVABasedHaloTagger.cc @@ -0,0 +1,468 @@ +/** \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 + */ + +#include "RecoEgamma/PhotonIdentification/interface/PhotonMVABasedHaloTagger.h" +#include "DataFormats/DetId/interface/DetId.h" +#include "DataFormats/EcalDetId/interface/EBDetId.h" +#include "DataFormats/EcalRecHit/interface/EcalRecHit.h" + +#include "FWCore/MessageLogger/interface/MessageLogger.h" +#include "RecoEcal/EgammaCoreTools/interface/EcalClusterTools.h" +#include +#include "CommonTools/MVAUtils/interface/GBRForestTools.h" + + +PhotonMVABasedHaloTagger::PhotonMVABasedHaloTagger(const edm::ParameterSet& conf, edm::ConsumesCollector&& iC) + : geometryToken_(iC.esConsumes()), + ecalPFRechitThresholdsToken_(iC.esConsumes()), + ecalClusterToolsESGetTokens_(std::move(iC)) +{ + trainingFileName_ = conf.getParameter("trainingFileName"); + rhoLabel_ = iC.consumes(conf.getParameter("rhoLabel")); + + EBecalCollection_ = iC.consumes(conf.getParameter("barrelEcalRecHitCollection")); + EEecalCollection_ = iC.consumes(conf.getParameter("endcapEcalRecHitCollection")); + ESCollection_ = iC.consumes(conf.getParameter("esRecHitCollection")); + HBHERecHitsCollection_ = iC.consumes(conf.getParameter("HBHERecHitsCollection")); + recHitEThresholdHB_ = conf.getParameter("recHitEThresholdHB"); + recHitEThresholdHE_ = conf.getParameter("recHitEThresholdHE"); + + noiseThrES_ = conf.getParameter("noiseThrES"); + gbr_ = createGBRForest(trainingFileName_); + +} + +double PhotonMVABasedHaloTagger::calculateMVA(const reco::Photon* pho, const edm::Event& iEvent, const edm::EventSetup& es){ + + + bool isEB = pho->isEB(); + + if(isEB) return 1.0; /// this MVA is useful and trained only for the EE photons. For EB, there are a lot of other useful handles which can reject beam halo efficiently + + //rho handle + edm::Handle rhoHandle; + iEvent.getByToken(rhoLabel_, rhoHandle); + + double rho_ = *(rhoHandle.product()); + + // Get all the RecHits + edm::Handle barrelHitHandle; + iEvent.getByToken(EBecalCollection_, barrelHitHandle); + + edm::Handle endcapHitHandle; + iEvent.getByToken(EEecalCollection_, endcapHitHandle); + + edm::Handle esHitHandle; + iEvent.getByToken(ESCollection_, esHitHandle); + + edm::Handle hbheHitHandle; + iEvent.getByToken(HBHERecHitsCollection_, hbheHitHandle); + + //gets geometry + pG_ = es.getHandle(geometryToken_); + const CaloGeometry* geo = pG_.product(); + + ///ECAL PF rechit thresholds + auto const& thresholds = es.getData(ecalPFRechitThresholdsToken_); + + noZS::EcalClusterLazyTools lazyToolnoZS(iEvent, ecalClusterToolsESGetTokens_.get(es), EBecalCollection_, EEecalCollection_); + + + ///calculate the energy weighted X, Y and Z position of the photon cluster + if(isEB) calphoClusCoordinECAL(geo, pho, &thresholds, barrelHitHandle); + else calphoClusCoordinECAL(geo, pho, &thresholds, endcapHitHandle); + + ///calculate the HBHE cluster position hypothesis + calmatchedHBHECoordForBothHypothesis(geo, pho, hbheHitHandle); + calmatchedESCoordForBothHypothesis(geo, pho, esHitHandle); + + ///this function works for EE only. Above ones work for EB as well in case later one wants to put a similar function for EB without returning 1 + + double angle_EE_HE_samedPhi = calAngleBetweenEEAndSubDet(hcalClusNhits_samedPhi_, hcalClusX_samedPhi_, hcalClusY_samedPhi_, hcalClusZ_samedPhi_); //essentially caculates the angle and energy variables in the two hypothesis between EE and HE + + double angle_EE_HE_samedR = calAngleBetweenEEAndSubDet(hcalClusNhits_samedR_, hcalClusX_samedR_, hcalClusY_samedR_, hcalClusZ_samedR_); + + double angle_EE_ES_samedPhi = calAngleBetweenEEAndSubDet(preshowerNhits_samedPhi_, preshowerX_samedPhi_, preshowerY_samedPhi_, preshowerZ_samedPhi_); + + double angle_EE_ES_samedR = calAngleBetweenEEAndSubDet(preshowerNhits_samedR_, preshowerX_samedR_, preshowerY_samedR_, preshowerZ_samedR_); + + ////set all the above calculated variables as input to the MVA + + const auto& vCov = lazyToolnoZS.localCovariances( *(pho->superCluster()->seed()) ); + double spp = (isnan(vCov[2]) ? 0. : sqrt(vCov[2])); + + ///https://cmssdt.cern.ch/lxr/source/RecoEgamma/ElectronIdentification/src/ElectronMVAEstimator.cc + + float vars[15]; + + vars[0] = preshowerE_samedPhi_; + vars[1] = hcalClusE_samedPhi_; + vars[2] = preshowerE_samedR_; + vars[3] = hcalClusE_samedR_; + vars[4] = pho->full5x5_r9(); + vars[5] = pho->superCluster()->etaWidth(); + vars[6] = pho->superCluster()->phiWidth(); + vars[7] = pho->full5x5_sigmaIetaIeta(); + vars[8] = spp; + vars[9] = angle_EE_ES_samedR; + vars[10] = angle_EE_HE_samedR; + vars[11] = angle_EE_ES_samedPhi; + vars[12] = angle_EE_HE_samedPhi; + vars[13] = (pho->superCluster()->preshowerEnergyPlane1() + pho->superCluster()->preshowerEnergyPlane2())/pho->superCluster()->rawEnergy(); + vars[14] = rho_; + + double BHmva = gbr_->GetAdaBoostClassifier(vars); + return BHmva; + +} + + +void PhotonMVABasedHaloTagger::calphoClusCoordinECAL(const CaloGeometry* geo, + const reco::Photon* pho, + const EcalPFRecHitThresholds *thresholds, + edm::Handle ecalrechitCollHandle + ){ + + ecalClusX_ = 0; ecalClusY_ = 0; ecalClusZ_ = 0; + ecalClusNhits_ = 0; + ecalClusE_ = 0; + + double phoSCEta = pho->superCluster()->eta(); + double phoSCPhi = pho->superCluster()->phi(); + + const EcalRecHitCollection* ecalRecHits = ecalrechitCollHandle.product(); + + for(EcalRecHitCollection::const_iterator ecalrechit = ecalRecHits->begin(); ecalrechit != ecalRecHits->end(); ecalrechit++){ + + auto const det = ecalrechit->id(); + double rhE = ecalrechit->energy(); + const GlobalPoint &rechitPoint = geo->getPosition(det); + + double rhEta = rechitPoint.eta(); + double rhPhi = rechitPoint.phi(); + double rhX = rechitPoint.x(); + double rhY = rechitPoint.y(); + double rhZ = rechitPoint.z(); + + if ((thresholds == nullptr)) { + throw cms::Exception("EmptyPFRechHitThresCollection") + << "In PhotonMVABasedHaloTagger::calphoClusCoordinECAL, EcalPFRecHitThresholds cannot be = nulptr"; + } + + float rhThres = 0.0; + if (thresholds != nullptr) { + rhThres = (*thresholds)[det]; + } + + if (rhE <= rhThres) + continue; + + + if(phoSCEta * rhEta < 0) continue; + + //double rho = sqrt(pow(rhX,2) + pow(rhY,2)); + double dPhi = deltaPhi(rhPhi, phoSCPhi); + double dEta = fabs(rhEta - phoSCEta); + + double dR = sqrt( pow(dPhi,2) + pow(dEta,2) ); + + + if(dR < 0.2){ + + ecalClusX_ += rhX * rhE; + ecalClusY_ += rhY * rhE; + ecalClusZ_ += rhZ * rhE; + ecalClusE_ += rhE; + ecalClusNhits_++; + + } + }//for(int ih=0; ih 0){ //should always be > 0 for an EM cluster + ecalClusX_ = ecalClusX_/ecalClusE_; + ecalClusY_ = ecalClusY_/ecalClusE_; + ecalClusZ_ = ecalClusZ_/ecalClusE_; + }//if(ecalClusNhits_>0) + + +} + +void PhotonMVABasedHaloTagger::calmatchedHBHECoordForBothHypothesis( + const CaloGeometry* geo, + const reco::Photon* pho, + edm::Handle hcalhitsCollHBHE + ){ + + + hcalClusX_samedPhi_ = 0; hcalClusY_samedPhi_ = 0; hcalClusZ_samedPhi_ = 0; + + hcalClusNhits_samedPhi_ = 0; + hcalClusE_samedPhi_ = 0; + + hcalClusX_samedR_ = 0; hcalClusY_samedR_ = 0; hcalClusZ_samedR_ = 0; + hcalClusNhits_samedR_ = 0; + hcalClusE_samedR_ = 0; + + const HBHERecHitCollection* HBHERecHits = hcalhitsCollHBHE.product(); + + double phoSCEta = pho->superCluster()->eta(); + double phoSCPhi = pho->superCluster()->phi(); + + + // Loop over HBHERecHit's + for (HBHERecHitCollection::const_iterator hbherechit = HBHERecHits->begin(); hbherechit != HBHERecHits->end(); hbherechit++) { + + HcalDetId det = hbherechit->id(); + const GlobalPoint & rechitPoint = geo->getPosition(det); + + double rhEta = rechitPoint.eta(); + double rhPhi = rechitPoint.phi(); + double rhX = rechitPoint.x(); + double rhY = rechitPoint.y(); + double rhZ = rechitPoint.z(); + double rhE = hbherechit->energy(); + + int depth = det.depth(); + + if ((det.subdet() == HcalBarrel and (depth < 1 or depth > int(recHitEThresholdHB_.size()))) or + (det.subdet() == HcalEndcap and (depth < 1 or depth > int(recHitEThresholdHE_.size())))) + edm::LogWarning("PhotonMVABasedHaloTagger") + << " hit in subdet " << det.subdet() << " has an unaccounted for depth of " << depth << "!!"; + + const bool goodHBe = det.subdet() == HcalBarrel and rhE > recHitEThresholdHB_[depth-1]; + const bool goodHEe = det.subdet() == HcalEndcap and rhE > recHitEThresholdHE_[depth-1]; + if (!(goodHBe or goodHEe)) + continue; + + if(phoSCEta * rhEta < 0) continue; ///Should be on the same side of Z + + + + double rho = sqrt(pow(rhX,2) + pow(rhY,2)); + double dEta = fabs(phoSCEta - rhEta); + double dPhi = deltaPhi(phoSCPhi, rhPhi); + + + bool isRHBehindECAL = false; + + double dRho = sqrt( pow(rhX-ecalClusX_,2) + pow(rhY-ecalClusY_,2) ) ; + + + if(rho>=31 && rho<=172 && dRho <= 26 && fabs(dPhi)<0.15){ ///only valid for the EE; this is 26 cm; hit within 3x3 of HCAL centered at the EECAL xtal + hcalClusX_samedPhi_ += rhX * rhE; + hcalClusY_samedPhi_ += rhY * rhE; + hcalClusZ_samedPhi_ += rhZ * rhE; + hcalClusE_samedPhi_ += rhE; + hcalClusNhits_samedPhi_++; + isRHBehindECAL = true; + + }//if(rho>=31 && rho<=172) + + double dR = sqrt( pow(dEta,2) + pow(dPhi,2) ); + + if(dR<0.15 && !isRHBehindECAL){ ///dont use hits which are just behind the ECAL in the same phi region + hcalClusX_samedR_ += rhX * rhE; + hcalClusY_samedR_ += rhY * rhE; + hcalClusZ_samedR_ += rhZ * rhE; + hcalClusE_samedR_ += rhE; + hcalClusNhits_samedR_++; + + } + + + }//for(int ih=0; ih0){ + hcalClusX_samedPhi_ = hcalClusX_samedPhi_/hcalClusE_samedPhi_; + hcalClusY_samedPhi_ = hcalClusY_samedPhi_/hcalClusE_samedPhi_; + hcalClusZ_samedPhi_ = hcalClusZ_samedPhi_/hcalClusE_samedPhi_; + }//if(hcalClusNhits_samedPhi_>0) + + if(hcalClusNhits_samedR_>0){ + hcalClusX_samedR_ = hcalClusX_samedR_/hcalClusE_samedR_; + hcalClusY_samedR_ = hcalClusY_samedR_/hcalClusE_samedR_; + hcalClusZ_samedR_ = hcalClusZ_samedR_/hcalClusE_samedR_; + }//if(hcalClusNhits_samedR_>0) +} + +void PhotonMVABasedHaloTagger::calmatchedESCoordForBothHypothesis(const CaloGeometry* geo, + const reco::Photon *pho, + edm::Handle esrechitCollHandle + ){ + + preshowerX_samedPhi_ = 0; preshowerY_samedPhi_ = 0; preshowerZ_samedPhi_ = 0; + preshowerNhits_samedPhi_ = 0; + preshowerE_samedPhi_ = 0; + + preshowerX_samedR_ = 0; preshowerY_samedR_ = 0; preshowerZ_samedR_ = 0; + preshowerNhits_samedR_ = 0; + preshowerE_samedR_ = 0; + + const EcalRecHitCollection* ESRecHits = esrechitCollHandle.product(); + + double phoSCEta = pho->superCluster()->eta(); + double phoSCPhi = pho->superCluster()->phi(); + + double tmpDiffdRho = 999; + double matchX_samephi = -999; + double matchY_samephi = -999; + bool foundESRH_samephi = false; + + double tmpDiffdRho_samedR = 999; + double matchX_samedR = -999; + double matchY_samedR = -999; + bool foundESRH_samedR = false; + + ///get theta and phi of the coordinates of photon + double tan_theta = 1./sinh(phoSCEta); + + double cos_phi = cos(phoSCPhi); + double sin_phi = sin(phoSCPhi); + + for( EcalRecHitCollection::const_iterator esrechit = ESRecHits->begin(); esrechit != ESRecHits->end(); esrechit++ ){ + + + const GlobalPoint & rechitPoint = geo->getPosition(esrechit->id()); + + double rhEta = rechitPoint.eta(); + //double rhPhi = rechitPoint.phi(); + double rhX = rechitPoint.x(); + double rhY = rechitPoint.y(); + double rhZ = rechitPoint.z(); + double rhE = esrechit->energy(); + + + if(phoSCEta * rhEta < 0) continue; + + + if(rhE < noiseThrES_) continue; + + /*double rho = sqrt(pow(rhX,2) + pow(rhY,2)); + double dEta = fabs(phoSCEta - rhEta); + double dPhi = deltaPhi(phoSCPhi, rhPhi); + */ + + ////try to include RH according to the strips, 11 in X and 11 in Y + /////////First calculate RH nearest in phi and eta to that of the photon SC + + //////same phi ----> the X and Y should be similar + double dRho = sqrt( pow(rhX-ecalClusX_,2) + pow(rhY-ecalClusY_,2) ) ; + + if(dRho < tmpDiffdRho && dRho < 2.2){ ////i.e. hit is required to be within the ----> seems better match with the data compared to 2.47 + + tmpDiffdRho = dRho; + matchX_samephi = rhX; + matchY_samephi = rhY; + foundESRH_samephi = true; + } + + + ////////same eta + ///calculate the expected x and y at the position of hte rechit + double exp_ESRho = rhZ * tan_theta; + double exp_ESX = cos_phi * exp_ESRho; + double exp_ESY = sin_phi * exp_ESRho; + + double dRho_samedR = sqrt( pow(rhX-exp_ESX, 2) + pow(rhY-exp_ESY,2) ); + + + if(dRho_samedR < tmpDiffdRho_samedR){ + + tmpDiffdRho_samedR = dRho_samedR; + matchX_samedR = rhX; + matchY_samedR = rhY; + foundESRH_samedR = true; + } + + }/// for( EcalRecHitCollection::const_iterator esrechit = ESRecHits->begin(); esrechit != ESRecHits->end(); esrechit++ ){ + + ////Now calculate the sum in +/- 5 strips in X and y around the matched RH + //+/5 strips mean = 5*~2mm = +/-10 mm = 1 cm + + for( EcalRecHitCollection::const_iterator esrechit = ESRecHits->begin(); esrechit != ESRecHits->end(); esrechit++ ){ + + const GlobalPoint & rechitPoint = geo->getPosition(esrechit->id()); + + double rhEta = rechitPoint.eta(); + double rhPhi = rechitPoint.phi(); + double rhX = rechitPoint.x(); + double rhY = rechitPoint.y(); + double rhZ = rechitPoint.z(); + double rhE = esrechit->energy(); + + + if(phoSCEta * rhEta < 0) continue; + if(rhE < noiseThrES_) continue; + + + bool isRHBehindECAL = false; + + ///same phi + double dX_samephi = fabs(matchX_samephi - rhX); + double dY_samephi = fabs(matchY_samephi - rhY); + + if( (dX_samephi < 1 && dY_samephi < 1) && foundESRH_samephi){ + + preshowerX_samedPhi_ += rhX * rhE; + preshowerY_samedPhi_ += rhY * rhE; + preshowerZ_samedPhi_ += rhZ * rhE; + preshowerE_samedPhi_ += rhE; + preshowerNhits_samedPhi_++; + isRHBehindECAL = true; + + } + + ///same dR + double dX_samedR = fabs(matchX_samedR - rhX); + double dY_samedR = fabs(matchY_samedR - rhY); + + if(!isRHBehindECAL && foundESRH_samedR && (dX_samedR < 1 && dY_samedR < 1) ){ + preshowerX_samedR_ += rhX * rhE; + preshowerY_samedR_ += rhY * rhE; + preshowerZ_samedR_ += rhZ * rhE; + preshowerE_samedR_ += rhE; + preshowerNhits_samedR_++; + } + + + }///for(int ih=0; ih0){ + preshowerX_samedPhi_ = preshowerX_samedPhi_/preshowerE_samedPhi_; + preshowerY_samedPhi_ = preshowerY_samedPhi_/preshowerE_samedPhi_; + preshowerZ_samedPhi_ = preshowerZ_samedPhi_/preshowerE_samedPhi_; + }//if(preshowerNhits_samedPhi_>0) + + + if(preshowerNhits_samedR_>0){ + preshowerX_samedR_ = preshowerX_samedR_/preshowerE_samedR_; + preshowerY_samedR_ = preshowerY_samedR_/preshowerE_samedR_; + preshowerZ_samedR_ = preshowerZ_samedR_/preshowerE_samedR_; + }//if(preshowerNhits_samedR_>0) +} + + +double PhotonMVABasedHaloTagger::calAngleBetweenEEAndSubDet(int nhits, double subdetClusX, double subdetClusY, double subdetClusZ){ + + ////get the angle of the line joining the ECAL cluster and the subdetector wrt Z axis for any hypothesis + + double angle = -999; + + if(ecalClusNhits_>0 && nhits>0 ){ + + double dR = sqrt( pow(subdetClusX-ecalClusX_,2) + pow(subdetClusY-ecalClusY_,2) + pow(subdetClusZ-ecalClusZ_,2) ); + + double cosTheta = fabs(subdetClusZ-ecalClusZ_)/dR; + + angle = acos(cosTheta); + } + + return angle; + +} + + + From ee342ca118b2ea6580b21575a5c2711b2ab7df04 Mon Sep 17 00:00:00 2001 From: Shilpi Date: Wed, 2 Feb 2022 14:47:13 +0100 Subject: [PATCH 02/13] added missing parameter for PhotonProducer --- RecoEgamma/EgammaPhotonProducers/python/photons_cfi.py | 1 + 1 file changed, 1 insertion(+) diff --git a/RecoEgamma/EgammaPhotonProducers/python/photons_cfi.py b/RecoEgamma/EgammaPhotonProducers/python/photons_cfi.py index bc05e54af71c8..99411cf896367 100644 --- a/RecoEgamma/EgammaPhotonProducers/python/photons_cfi.py +++ b/RecoEgamma/EgammaPhotonProducers/python/photons_cfi.py @@ -164,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), From fdcfde27859c5f281ce8be5507d21536db444b5d Mon Sep 17 00:00:00 2001 From: Shilpi Date: Mon, 7 Feb 2022 01:27:35 +0100 Subject: [PATCH 03/13] Code format checks --- .../EgammaCandidates/interface/Photon.h | 6 +- DataFormats/EgammaCandidates/src/Photon.cc | 10 +- .../src/GEDPhotonProducer.cc | 6 +- .../src/PhotonProducer.cc | 6 +- .../interface/PhotonMVABasedHaloTagger.h | 58 +-- .../src/PhotonMVABasedHaloTagger.cc | 402 +++++++++--------- 6 files changed, 231 insertions(+), 257 deletions(-) diff --git a/DataFormats/EgammaCandidates/interface/Photon.h b/DataFormats/EgammaCandidates/interface/Photon.h index f4b3bc4e2bc1c..ddae30524c9b7 100644 --- a/DataFormats/EgammaCandidates/interface/Photon.h +++ b/DataFormats/EgammaCandidates/interface/Photon.h @@ -581,11 +581,9 @@ namespace reco { ///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; - } + void setHaloTaggerMVAVal(float x) { haloTaggerMVAVal_ = x; } private: /// check overlap with another candidate diff --git a/DataFormats/EgammaCandidates/src/Photon.cc b/DataFormats/EgammaCandidates/src/Photon.cc index 27892904b9144..7633e35a966f5 100644 --- a/DataFormats/EgammaCandidates/src/Photon.cc +++ b/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),haloTaggerMVAVal_(-999) {} + : RecoCandidate(0, p4, vtx, 22), + caloPosition_(caloPos), + photonCore_(core), + pixelSeed_(false), + haloTaggerMVAVal_(-999) {} Photon::Photon(const Photon& rhs) : RecoCandidate(rhs), @@ -19,8 +23,8 @@ Photon::Photon(const Photon& rhs) saturationInfo_(rhs.saturationInfo_), eCorrections_(rhs.eCorrections_), mipVariableBlock_(rhs.mipVariableBlock_), - pfIsolation_(rhs.pfIsolation_), - haloTaggerMVAVal_(rhs.haloTaggerMVAVal_) {} + pfIsolation_(rhs.pfIsolation_), + haloTaggerMVAVal_(rhs.haloTaggerMVAVal_) {} Photon::~Photon() {} diff --git a/RecoEgamma/EgammaPhotonProducers/src/GEDPhotonProducer.cc b/RecoEgamma/EgammaPhotonProducers/src/GEDPhotonProducer.cc index 30713677ace12..b85643ac35874 100644 --- a/RecoEgamma/EgammaPhotonProducers/src/GEDPhotonProducer.cc +++ b/RecoEgamma/EgammaPhotonProducers/src/GEDPhotonProducer.cc @@ -418,8 +418,8 @@ GEDPhotonProducer::GEDPhotonProducer(const edm::ParameterSet& config, const Cach if (recoStep_.isFinal() && runMVABasedHaloTagger_) { edm::ParameterSet mvaBasedHaloVariableSet = config.getParameter("mvaBasedHaloVariableSet"); - photonMVABasedHaloTagger_ = std::make_unique(mvaBasedHaloVariableSet, consumesCollector()); - + photonMVABasedHaloTagger_ = + std::make_unique(mvaBasedHaloVariableSet, consumesCollector()); } ///Get the set for PF cluster isolation calculator @@ -1068,7 +1068,7 @@ void GEDPhotonProducer::fillPhotonCollection(edm::Event& evt, reco::Photon newCandidate(*phoRef); iSC++; - if (runMVABasedHaloTagger_) { ///sets values only for EE, for EB it always returns 1 + if (runMVABasedHaloTagger_) { ///sets values only for EE, for EB it always returns 1 float BHmva = photonMVABasedHaloTagger_->calculateMVA(&newCandidate, evt, es); newCandidate.setHaloTaggerMVAVal(BHmva); } diff --git a/RecoEgamma/EgammaPhotonProducers/src/PhotonProducer.cc b/RecoEgamma/EgammaPhotonProducers/src/PhotonProducer.cc index 4db6d3e9b0074..e733bc925bad5 100644 --- a/RecoEgamma/EgammaPhotonProducers/src/PhotonProducer.cc +++ b/RecoEgamma/EgammaPhotonProducers/src/PhotonProducer.cc @@ -238,9 +238,9 @@ PhotonProducer::PhotonProducer(const edm::ParameterSet& config) if (runMVABasedHaloTagger_) { edm::ParameterSet mvaBasedHaloVariableSet = config.getParameter("mvaBasedHaloVariableSet"); - photonMVABasedHaloTagger_ = std::make_unique(mvaBasedHaloVariableSet, consumesCollector()); + photonMVABasedHaloTagger_ = + std::make_unique(mvaBasedHaloVariableSet, consumesCollector()); } - // Register the product produces(PhotonCollection_); @@ -498,7 +498,7 @@ void PhotonProducer::fillPhotonCollection(edm::Event& evt, newCandidate.setMIPVariables(mipVar); } - if (runMVABasedHaloTagger_) { ///sets values only for EE, for EB it always returns 1 + if (runMVABasedHaloTagger_) { ///sets values only for EE, for EB it always returns 1 double BHmva = photonMVABasedHaloTagger_->calculateMVA(&newCandidate, evt, es); newCandidate.setHaloTaggerMVAVal(BHmva); } diff --git a/RecoEgamma/PhotonIdentification/interface/PhotonMVABasedHaloTagger.h b/RecoEgamma/PhotonIdentification/interface/PhotonMVABasedHaloTagger.h index f3498f526fc8b..b25d0e2f3e35f 100644 --- a/RecoEgamma/PhotonIdentification/interface/PhotonMVABasedHaloTagger.h +++ b/RecoEgamma/PhotonIdentification/interface/PhotonMVABasedHaloTagger.h @@ -20,7 +20,7 @@ #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/EcalObjects/interface/EcalPFRecHitThresholds.h" #include "CondFormats/DataRecord/interface/EcalPFRecHitThresholdsRcd.h" #include @@ -29,64 +29,52 @@ class PhotonMVABasedHaloTagger { PhotonMVABasedHaloTagger(const edm::ParameterSet& conf, edm::ConsumesCollector&& iC); virtual ~PhotonMVABasedHaloTagger(){}; - double calculateMVA(const reco::Photon* pho, - const edm::Event&iEvent, - const edm::EventSetup& es - ); - - + double calculateMVA(const reco::Photon* pho, const edm::Event& iEvent, const edm::EventSetup& es); + void calphoClusCoordinECAL(const CaloGeometry* geo, - const reco::Photon*, - const EcalPFRecHitThresholds *thresholds, - edm::Handle ecalrechitCollHandle - ); - + const reco::Photon*, + const EcalPFRecHitThresholds* thresholds, + edm::Handle ecalrechitCollHandle); void calmatchedHBHECoordForBothHypothesis(const CaloGeometry* geo, - const reco::Photon*, - edm::Handle hbheHitHandle - ); - - + const reco::Photon*, + edm::Handle hbheHitHandle); + void calmatchedESCoordForBothHypothesis(const CaloGeometry* geo, - const reco::Photon*, - edm::Handle esrechitCollHandle - ); - - double calAngleBetweenEEAndSubDet(int nhits, - double subdetClusX, - double subdetClusY, - double subdetClusZ); - + const reco::Photon*, + edm::Handle esrechitCollHandle); + + double calAngleBetweenEEAndSubDet(int nhits, double subdetClusX, double subdetClusY, double subdetClusZ); private: int hcalClusNhits_samedPhi_, hcalClusNhits_samedR_; - int ecalClusNhits_, preshowerNhits_samedPhi_, preshowerNhits_samedR_; - double hcalClusX_samedPhi_, hcalClusY_samedPhi_, hcalClusZ_samedPhi_, hcalClusX_samedR_, hcalClusY_samedR_, hcalClusZ_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 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_; - std::string trainingFileName_; + std::string trainingFileName_; const edm::ESGetToken geometryToken_; const edm::ESGetToken ecalPFRechitThresholdsToken_; const EcalClusterLazyTools::ESGetTokens ecalClusterToolsESGetTokens_; - - edm::ESHandle pG_; - edm::EDGetTokenT rhoLabel_; + + edm::ESHandle pG_; + edm::EDGetTokenT rhoLabel_; edm::EDGetTokenT EBecalCollection_; edm::EDGetTokenT EEecalCollection_; edm::EDGetTokenT ESCollection_; edm::EDGetTokenT HBHERecHitsCollection_; std::unique_ptr gbr_; - }; #endif // PhotonMVABasedHaloTagger_H diff --git a/RecoEgamma/PhotonIdentification/src/PhotonMVABasedHaloTagger.cc b/RecoEgamma/PhotonIdentification/src/PhotonMVABasedHaloTagger.cc index 9154c8a7efc9c..29fdc57e3c23f 100644 --- a/RecoEgamma/PhotonIdentification/src/PhotonMVABasedHaloTagger.cc +++ b/RecoEgamma/PhotonIdentification/src/PhotonMVABasedHaloTagger.cc @@ -15,15 +15,13 @@ #include #include "CommonTools/MVAUtils/interface/GBRForestTools.h" - PhotonMVABasedHaloTagger::PhotonMVABasedHaloTagger(const edm::ParameterSet& conf, edm::ConsumesCollector&& iC) - : geometryToken_(iC.esConsumes()), - ecalPFRechitThresholdsToken_(iC.esConsumes()), - ecalClusterToolsESGetTokens_(std::move(iC)) -{ + : geometryToken_(iC.esConsumes()), + ecalPFRechitThresholdsToken_(iC.esConsumes()), + ecalClusterToolsESGetTokens_(std::move(iC)) { trainingFileName_ = conf.getParameter("trainingFileName"); - rhoLabel_ = iC.consumes(conf.getParameter("rhoLabel")); - + rhoLabel_ = iC.consumes(conf.getParameter("rhoLabel")); + EBecalCollection_ = iC.consumes(conf.getParameter("barrelEcalRecHitCollection")); EEecalCollection_ = iC.consumes(conf.getParameter("endcapEcalRecHitCollection")); ESCollection_ = iC.consumes(conf.getParameter("esRecHitCollection")); @@ -31,34 +29,34 @@ PhotonMVABasedHaloTagger::PhotonMVABasedHaloTagger(const edm::ParameterSet& conf recHitEThresholdHB_ = conf.getParameter("recHitEThresholdHB"); recHitEThresholdHE_ = conf.getParameter("recHitEThresholdHE"); - noiseThrES_ = conf.getParameter("noiseThrES"); + noiseThrES_ = conf.getParameter("noiseThrES"); gbr_ = createGBRForest(trainingFileName_); - } -double PhotonMVABasedHaloTagger::calculateMVA(const reco::Photon* pho, const edm::Event& iEvent, const edm::EventSetup& es){ - - +double PhotonMVABasedHaloTagger::calculateMVA(const reco::Photon* pho, + const edm::Event& iEvent, + const edm::EventSetup& es) { bool isEB = pho->isEB(); - if(isEB) return 1.0; /// this MVA is useful and trained only for the EE photons. For EB, there are a lot of other useful handles which can reject beam halo efficiently - + if (isEB) + return 1.0; /// this MVA is useful and trained only for the EE photons. For EB, there are a lot of other useful handles which can reject beam halo efficiently + //rho handle edm::Handle rhoHandle; iEvent.getByToken(rhoLabel_, rhoHandle); - double rho_ = *(rhoHandle.product()); + double rho_ = *(rhoHandle.product()); // Get all the RecHits edm::Handle barrelHitHandle; iEvent.getByToken(EBecalCollection_, barrelHitHandle); - + edm::Handle endcapHitHandle; iEvent.getByToken(EEecalCollection_, endcapHitHandle); edm::Handle esHitHandle; iEvent.getByToken(ESCollection_, esHitHandle); - + edm::Handle hbheHitHandle; iEvent.getByToken(HBHERecHitsCollection_, hbheHitHandle); @@ -68,33 +66,42 @@ double PhotonMVABasedHaloTagger::calculateMVA(const reco::Photon* pho, const edm ///ECAL PF rechit thresholds auto const& thresholds = es.getData(ecalPFRechitThresholdsToken_); - - noZS::EcalClusterLazyTools lazyToolnoZS(iEvent, ecalClusterToolsESGetTokens_.get(es), EBecalCollection_, EEecalCollection_); + noZS::EcalClusterLazyTools lazyToolnoZS( + iEvent, ecalClusterToolsESGetTokens_.get(es), EBecalCollection_, EEecalCollection_); ///calculate the energy weighted X, Y and Z position of the photon cluster - if(isEB) calphoClusCoordinECAL(geo, pho, &thresholds, barrelHitHandle); - else calphoClusCoordinECAL(geo, pho, &thresholds, endcapHitHandle); - + if (isEB) + calphoClusCoordinECAL(geo, pho, &thresholds, barrelHitHandle); + else + calphoClusCoordinECAL(geo, pho, &thresholds, endcapHitHandle); + ///calculate the HBHE cluster position hypothesis calmatchedHBHECoordForBothHypothesis(geo, pho, hbheHitHandle); calmatchedESCoordForBothHypothesis(geo, pho, esHitHandle); ///this function works for EE only. Above ones work for EB as well in case later one wants to put a similar function for EB without returning 1 - double angle_EE_HE_samedPhi = calAngleBetweenEEAndSubDet(hcalClusNhits_samedPhi_, hcalClusX_samedPhi_, hcalClusY_samedPhi_, hcalClusZ_samedPhi_); //essentially caculates the angle and energy variables in the two hypothesis between EE and HE + double angle_EE_HE_samedPhi = calAngleBetweenEEAndSubDet( + hcalClusNhits_samedPhi_, + hcalClusX_samedPhi_, + hcalClusY_samedPhi_, + hcalClusZ_samedPhi_); //essentially caculates the angle and energy variables in the two hypothesis between EE and HE - double angle_EE_HE_samedR = calAngleBetweenEEAndSubDet(hcalClusNhits_samedR_, hcalClusX_samedR_, hcalClusY_samedR_, hcalClusZ_samedR_); + double angle_EE_HE_samedR = + calAngleBetweenEEAndSubDet(hcalClusNhits_samedR_, hcalClusX_samedR_, hcalClusY_samedR_, hcalClusZ_samedR_); - double angle_EE_ES_samedPhi = calAngleBetweenEEAndSubDet(preshowerNhits_samedPhi_, preshowerX_samedPhi_, preshowerY_samedPhi_, preshowerZ_samedPhi_); + double angle_EE_ES_samedPhi = calAngleBetweenEEAndSubDet( + preshowerNhits_samedPhi_, preshowerX_samedPhi_, preshowerY_samedPhi_, preshowerZ_samedPhi_); + + double angle_EE_ES_samedR = + calAngleBetweenEEAndSubDet(preshowerNhits_samedR_, preshowerX_samedR_, preshowerY_samedR_, preshowerZ_samedR_); - double angle_EE_ES_samedR = calAngleBetweenEEAndSubDet(preshowerNhits_samedR_, preshowerX_samedR_, preshowerY_samedR_, preshowerZ_samedR_); - ////set all the above calculated variables as input to the MVA - const auto& vCov = lazyToolnoZS.localCovariances( *(pho->superCluster()->seed()) ); + const auto& vCov = lazyToolnoZS.localCovariances(*(pho->superCluster()->seed())); double spp = (isnan(vCov[2]) ? 0. : sqrt(vCov[2])); - + ///https://cmssdt.cern.ch/lxr/source/RecoEgamma/ElectronIdentification/src/ElectronMVAEstimator.cc float vars[15]; @@ -112,221 +119,208 @@ double PhotonMVABasedHaloTagger::calculateMVA(const reco::Photon* pho, const edm vars[10] = angle_EE_HE_samedR; vars[11] = angle_EE_ES_samedPhi; vars[12] = angle_EE_HE_samedPhi; - vars[13] = (pho->superCluster()->preshowerEnergyPlane1() + pho->superCluster()->preshowerEnergyPlane2())/pho->superCluster()->rawEnergy(); + vars[13] = (pho->superCluster()->preshowerEnergyPlane1() + pho->superCluster()->preshowerEnergyPlane2()) / + pho->superCluster()->rawEnergy(); vars[14] = rho_; double BHmva = gbr_->GetAdaBoostClassifier(vars); return BHmva; - } - void PhotonMVABasedHaloTagger::calphoClusCoordinECAL(const CaloGeometry* geo, - const reco::Photon* pho, - const EcalPFRecHitThresholds *thresholds, - edm::Handle ecalrechitCollHandle - ){ - - ecalClusX_ = 0; ecalClusY_ = 0; ecalClusZ_ = 0; + const reco::Photon* pho, + const EcalPFRecHitThresholds* thresholds, + edm::Handle ecalrechitCollHandle) { + ecalClusX_ = 0; + ecalClusY_ = 0; + ecalClusZ_ = 0; ecalClusNhits_ = 0; ecalClusE_ = 0; double phoSCEta = pho->superCluster()->eta(); double phoSCPhi = pho->superCluster()->phi(); - + const EcalRecHitCollection* ecalRecHits = ecalrechitCollHandle.product(); - - for(EcalRecHitCollection::const_iterator ecalrechit = ecalRecHits->begin(); ecalrechit != ecalRecHits->end(); ecalrechit++){ - + + for (EcalRecHitCollection::const_iterator ecalrechit = ecalRecHits->begin(); ecalrechit != ecalRecHits->end(); + ecalrechit++) { auto const det = ecalrechit->id(); double rhE = ecalrechit->energy(); - const GlobalPoint &rechitPoint = geo->getPosition(det); - + const GlobalPoint& rechitPoint = geo->getPosition(det); + double rhEta = rechitPoint.eta(); double rhPhi = rechitPoint.phi(); double rhX = rechitPoint.x(); double rhY = rechitPoint.y(); double rhZ = rechitPoint.z(); - + if ((thresholds == nullptr)) { throw cms::Exception("EmptyPFRechHitThresCollection") - << "In PhotonMVABasedHaloTagger::calphoClusCoordinECAL, EcalPFRecHitThresholds cannot be = nulptr"; + << "In PhotonMVABasedHaloTagger::calphoClusCoordinECAL, EcalPFRecHitThresholds cannot be = nulptr"; } - + float rhThres = 0.0; if (thresholds != nullptr) { - rhThres = (*thresholds)[det]; + rhThres = (*thresholds)[det]; } - + if (rhE <= rhThres) continue; - - if(phoSCEta * rhEta < 0) continue; + if (phoSCEta * rhEta < 0) + continue; //double rho = sqrt(pow(rhX,2) + pow(rhY,2)); double dPhi = deltaPhi(rhPhi, phoSCPhi); double dEta = fabs(rhEta - phoSCEta); - - double dR = sqrt( pow(dPhi,2) + pow(dEta,2) ); - - if(dR < 0.2){ - + double dR = sqrt(pow(dPhi, 2) + pow(dEta, 2)); + + if (dR < 0.2) { ecalClusX_ += rhX * rhE; ecalClusY_ += rhY * rhE; ecalClusZ_ += rhZ * rhE; ecalClusE_ += rhE; ecalClusNhits_++; - } - }//for(int ih=0; ih 0){ //should always be > 0 for an EM cluster - ecalClusX_ = ecalClusX_/ecalClusE_; - ecalClusY_ = ecalClusY_/ecalClusE_; - ecalClusZ_ = ecalClusZ_/ecalClusE_; - }//if(ecalClusNhits_>0) - + } //for(int ih=0; ih 0) { //should always be > 0 for an EM cluster + ecalClusX_ = ecalClusX_ / ecalClusE_; + ecalClusY_ = ecalClusY_ / ecalClusE_; + ecalClusZ_ = ecalClusZ_ / ecalClusE_; + } //if(ecalClusNhits_>0) } void PhotonMVABasedHaloTagger::calmatchedHBHECoordForBothHypothesis( - const CaloGeometry* geo, - const reco::Photon* pho, - edm::Handle hcalhitsCollHBHE - ){ - - - hcalClusX_samedPhi_ = 0; hcalClusY_samedPhi_ = 0; hcalClusZ_samedPhi_ = 0; - + const CaloGeometry* geo, const reco::Photon* pho, edm::Handle hcalhitsCollHBHE) { + hcalClusX_samedPhi_ = 0; + hcalClusY_samedPhi_ = 0; + hcalClusZ_samedPhi_ = 0; + hcalClusNhits_samedPhi_ = 0; hcalClusE_samedPhi_ = 0; - - hcalClusX_samedR_ = 0; hcalClusY_samedR_ = 0; hcalClusZ_samedR_ = 0; + + hcalClusX_samedR_ = 0; + hcalClusY_samedR_ = 0; + hcalClusZ_samedR_ = 0; hcalClusNhits_samedR_ = 0; hcalClusE_samedR_ = 0; - - const HBHERecHitCollection* HBHERecHits = hcalhitsCollHBHE.product(); - + + const HBHERecHitCollection* HBHERecHits = hcalhitsCollHBHE.product(); + double phoSCEta = pho->superCluster()->eta(); double phoSCPhi = pho->superCluster()->phi(); - // Loop over HBHERecHit's - for (HBHERecHitCollection::const_iterator hbherechit = HBHERecHits->begin(); hbherechit != HBHERecHits->end(); hbherechit++) { - + for (HBHERecHitCollection::const_iterator hbherechit = HBHERecHits->begin(); hbherechit != HBHERecHits->end(); + hbherechit++) { HcalDetId det = hbherechit->id(); - const GlobalPoint & rechitPoint = geo->getPosition(det); - + const GlobalPoint& rechitPoint = geo->getPosition(det); + double rhEta = rechitPoint.eta(); double rhPhi = rechitPoint.phi(); double rhX = rechitPoint.x(); double rhY = rechitPoint.y(); double rhZ = rechitPoint.z(); double rhE = hbherechit->energy(); - + int depth = det.depth(); - + if ((det.subdet() == HcalBarrel and (depth < 1 or depth > int(recHitEThresholdHB_.size()))) or - (det.subdet() == HcalEndcap and (depth < 1 or depth > int(recHitEThresholdHE_.size())))) + (det.subdet() == HcalEndcap and (depth < 1 or depth > int(recHitEThresholdHE_.size())))) edm::LogWarning("PhotonMVABasedHaloTagger") - << " hit in subdet " << det.subdet() << " has an unaccounted for depth of " << depth << "!!"; - - const bool goodHBe = det.subdet() == HcalBarrel and rhE > recHitEThresholdHB_[depth-1]; - const bool goodHEe = det.subdet() == HcalEndcap and rhE > recHitEThresholdHE_[depth-1]; + << " hit in subdet " << det.subdet() << " has an unaccounted for depth of " << depth << "!!"; + + const bool goodHBe = det.subdet() == HcalBarrel and rhE > recHitEThresholdHB_[depth - 1]; + const bool goodHEe = det.subdet() == HcalEndcap and rhE > recHitEThresholdHE_[depth - 1]; if (!(goodHBe or goodHEe)) continue; - - if(phoSCEta * rhEta < 0) continue; ///Should be on the same side of Z - - - double rho = sqrt(pow(rhX,2) + pow(rhY,2)); + if (phoSCEta * rhEta < 0) + continue; ///Should be on the same side of Z + + double rho = sqrt(pow(rhX, 2) + pow(rhY, 2)); double dEta = fabs(phoSCEta - rhEta); double dPhi = deltaPhi(phoSCPhi, rhPhi); - - + bool isRHBehindECAL = false; - - double dRho = sqrt( pow(rhX-ecalClusX_,2) + pow(rhY-ecalClusY_,2) ) ; - - - if(rho>=31 && rho<=172 && dRho <= 26 && fabs(dPhi)<0.15){ ///only valid for the EE; this is 26 cm; hit within 3x3 of HCAL centered at the EECAL xtal + + double dRho = sqrt(pow(rhX - ecalClusX_, 2) + pow(rhY - ecalClusY_, 2)); + + if (rho >= 31 && rho <= 172 && dRho <= 26 && + fabs(dPhi) < 0.15) { ///only valid for the EE; this is 26 cm; hit within 3x3 of HCAL centered at the EECAL xtal hcalClusX_samedPhi_ += rhX * rhE; hcalClusY_samedPhi_ += rhY * rhE; hcalClusZ_samedPhi_ += rhZ * rhE; hcalClusE_samedPhi_ += rhE; hcalClusNhits_samedPhi_++; isRHBehindECAL = true; - - }//if(rho>=31 && rho<=172) - - double dR = sqrt( pow(dEta,2) + pow(dPhi,2) ); - - if(dR<0.15 && !isRHBehindECAL){ ///dont use hits which are just behind the ECAL in the same phi region + + } //if(rho>=31 && rho<=172) + + double dR = sqrt(pow(dEta, 2) + pow(dPhi, 2)); + + if (dR < 0.15 && !isRHBehindECAL) { ///dont use hits which are just behind the ECAL in the same phi region hcalClusX_samedR_ += rhX * rhE; hcalClusY_samedR_ += rhY * rhE; hcalClusZ_samedR_ += rhZ * rhE; hcalClusE_samedR_ += rhE; hcalClusNhits_samedR_++; - } - - - }//for(int ih=0; ih0){ - hcalClusX_samedPhi_ = hcalClusX_samedPhi_/hcalClusE_samedPhi_; - hcalClusY_samedPhi_ = hcalClusY_samedPhi_/hcalClusE_samedPhi_; - hcalClusZ_samedPhi_ = hcalClusZ_samedPhi_/hcalClusE_samedPhi_; - }//if(hcalClusNhits_samedPhi_>0) - - if(hcalClusNhits_samedR_>0){ - hcalClusX_samedR_ = hcalClusX_samedR_/hcalClusE_samedR_; - hcalClusY_samedR_ = hcalClusY_samedR_/hcalClusE_samedR_; - hcalClusZ_samedR_ = hcalClusZ_samedR_/hcalClusE_samedR_; - }//if(hcalClusNhits_samedR_>0) + + } //for(int ih=0; ih 0) { + hcalClusX_samedPhi_ = hcalClusX_samedPhi_ / hcalClusE_samedPhi_; + hcalClusY_samedPhi_ = hcalClusY_samedPhi_ / hcalClusE_samedPhi_; + hcalClusZ_samedPhi_ = hcalClusZ_samedPhi_ / hcalClusE_samedPhi_; + } //if(hcalClusNhits_samedPhi_>0) + + if (hcalClusNhits_samedR_ > 0) { + hcalClusX_samedR_ = hcalClusX_samedR_ / hcalClusE_samedR_; + hcalClusY_samedR_ = hcalClusY_samedR_ / hcalClusE_samedR_; + hcalClusZ_samedR_ = hcalClusZ_samedR_ / hcalClusE_samedR_; + } //if(hcalClusNhits_samedR_>0) } -void PhotonMVABasedHaloTagger::calmatchedESCoordForBothHypothesis(const CaloGeometry* geo, - const reco::Photon *pho, - edm::Handle esrechitCollHandle - ){ - - preshowerX_samedPhi_ = 0; preshowerY_samedPhi_ = 0; preshowerZ_samedPhi_ = 0; +void PhotonMVABasedHaloTagger::calmatchedESCoordForBothHypothesis( + const CaloGeometry* geo, const reco::Photon* pho, edm::Handle esrechitCollHandle) { + preshowerX_samedPhi_ = 0; + preshowerY_samedPhi_ = 0; + preshowerZ_samedPhi_ = 0; preshowerNhits_samedPhi_ = 0; preshowerE_samedPhi_ = 0; - preshowerX_samedR_ = 0; preshowerY_samedR_ = 0; preshowerZ_samedR_ = 0; + preshowerX_samedR_ = 0; + preshowerY_samedR_ = 0; + preshowerZ_samedR_ = 0; preshowerNhits_samedR_ = 0; preshowerE_samedR_ = 0; const EcalRecHitCollection* ESRecHits = esrechitCollHandle.product(); - + double phoSCEta = pho->superCluster()->eta(); double phoSCPhi = pho->superCluster()->phi(); - + double tmpDiffdRho = 999; double matchX_samephi = -999; double matchY_samephi = -999; bool foundESRH_samephi = false; - + double tmpDiffdRho_samedR = 999; double matchX_samedR = -999; double matchY_samedR = -999; bool foundESRH_samedR = false; ///get theta and phi of the coordinates of photon - double tan_theta = 1./sinh(phoSCEta); - + double tan_theta = 1. / sinh(phoSCEta); + double cos_phi = cos(phoSCPhi); double sin_phi = sin(phoSCPhi); - for( EcalRecHitCollection::const_iterator esrechit = ESRecHits->begin(); esrechit != ESRecHits->end(); esrechit++ ){ - - - const GlobalPoint & rechitPoint = geo->getPosition(esrechit->id()); + for (EcalRecHitCollection::const_iterator esrechit = ESRecHits->begin(); esrechit != ESRecHits->end(); esrechit++) { + const GlobalPoint& rechitPoint = geo->getPosition(esrechit->id()); double rhEta = rechitPoint.eta(); //double rhPhi = rechitPoint.phi(); @@ -334,58 +328,55 @@ void PhotonMVABasedHaloTagger::calmatchedESCoordForBothHypothesis(const CaloGeom double rhY = rechitPoint.y(); double rhZ = rechitPoint.z(); double rhE = esrechit->energy(); - - - if(phoSCEta * rhEta < 0) continue; + if (phoSCEta * rhEta < 0) + continue; + + if (rhE < noiseThrES_) + continue; - if(rhE < noiseThrES_) continue; - /*double rho = sqrt(pow(rhX,2) + pow(rhY,2)); double dEta = fabs(phoSCEta - rhEta); double dPhi = deltaPhi(phoSCPhi, rhPhi); */ - ////try to include RH according to the strips, 11 in X and 11 in Y - /////////First calculate RH nearest in phi and eta to that of the photon SC - + ////try to include RH according to the strips, 11 in X and 11 in Y + /////////First calculate RH nearest in phi and eta to that of the photon SC + //////same phi ----> the X and Y should be similar - double dRho = sqrt( pow(rhX-ecalClusX_,2) + pow(rhY-ecalClusY_,2) ) ; + double dRho = sqrt(pow(rhX - ecalClusX_, 2) + pow(rhY - ecalClusY_, 2)); + + if (dRho < tmpDiffdRho && + dRho < 2.2) { ////i.e. hit is required to be within the ----> seems better match with the data compared to 2.47 - if(dRho < tmpDiffdRho && dRho < 2.2){ ////i.e. hit is required to be within the ----> seems better match with the data compared to 2.47 - tmpDiffdRho = dRho; matchX_samephi = rhX; matchY_samephi = rhY; foundESRH_samephi = true; } - - ////////same eta - ///calculate the expected x and y at the position of hte rechit + ////////same eta + ///calculate the expected x and y at the position of hte rechit double exp_ESRho = rhZ * tan_theta; double exp_ESX = cos_phi * exp_ESRho; double exp_ESY = sin_phi * exp_ESRho; - double dRho_samedR = sqrt( pow(rhX-exp_ESX, 2) + pow(rhY-exp_ESY,2) ); - - - if(dRho_samedR < tmpDiffdRho_samedR){ + double dRho_samedR = sqrt(pow(rhX - exp_ESX, 2) + pow(rhY - exp_ESY, 2)); + if (dRho_samedR < tmpDiffdRho_samedR) { tmpDiffdRho_samedR = dRho_samedR; matchX_samedR = rhX; matchY_samedR = rhY; foundESRH_samedR = true; } - }/// for( EcalRecHitCollection::const_iterator esrechit = ESRecHits->begin(); esrechit != ESRecHits->end(); esrechit++ ){ - - ////Now calculate the sum in +/- 5 strips in X and y around the matched RH - //+/5 strips mean = 5*~2mm = +/-10 mm = 1 cm + } /// for( EcalRecHitCollection::const_iterator esrechit = ESRecHits->begin(); esrechit != ESRecHits->end(); esrechit++ ){ - for( EcalRecHitCollection::const_iterator esrechit = ESRecHits->begin(); esrechit != ESRecHits->end(); esrechit++ ){ + ////Now calculate the sum in +/- 5 strips in X and y around the matched RH + //+/5 strips mean = 5*~2mm = +/-10 mm = 1 cm - const GlobalPoint & rechitPoint = geo->getPosition(esrechit->id()); + for (EcalRecHitCollection::const_iterator esrechit = ESRecHits->begin(); esrechit != ESRecHits->end(); esrechit++) { + const GlobalPoint& rechitPoint = geo->getPosition(esrechit->id()); double rhEta = rechitPoint.eta(); double rhPhi = rechitPoint.phi(); @@ -393,76 +384,69 @@ void PhotonMVABasedHaloTagger::calmatchedESCoordForBothHypothesis(const CaloGeom double rhY = rechitPoint.y(); double rhZ = rechitPoint.z(); double rhE = esrechit->energy(); - - - if(phoSCEta * rhEta < 0) continue; - if(rhE < noiseThrES_) continue; + if (phoSCEta * rhEta < 0) + continue; + if (rhE < noiseThrES_) + continue; bool isRHBehindECAL = false; - + ///same phi double dX_samephi = fabs(matchX_samephi - rhX); double dY_samephi = fabs(matchY_samephi - rhY); - if( (dX_samephi < 1 && dY_samephi < 1) && foundESRH_samephi){ - + if ((dX_samephi < 1 && dY_samephi < 1) && foundESRH_samephi) { preshowerX_samedPhi_ += rhX * rhE; preshowerY_samedPhi_ += rhY * rhE; preshowerZ_samedPhi_ += rhZ * rhE; - preshowerE_samedPhi_ += rhE; + preshowerE_samedPhi_ += rhE; preshowerNhits_samedPhi_++; isRHBehindECAL = true; - } - ///same dR + ///same dR double dX_samedR = fabs(matchX_samedR - rhX); double dY_samedR = fabs(matchY_samedR - rhY); - if(!isRHBehindECAL && foundESRH_samedR && (dX_samedR < 1 && dY_samedR < 1) ){ - preshowerX_samedR_ += rhX * rhE; - preshowerY_samedR_ += rhY * rhE; - preshowerZ_samedR_ += rhZ * rhE; + if (!isRHBehindECAL && foundESRH_samedR && (dX_samedR < 1 && dY_samedR < 1)) { + preshowerX_samedR_ += rhX * rhE; + preshowerY_samedR_ += rhY * rhE; + preshowerZ_samedR_ += rhZ * rhE; preshowerE_samedR_ += rhE; preshowerNhits_samedR_++; } - - - }///for(int ih=0; ih0){ - preshowerX_samedPhi_ = preshowerX_samedPhi_/preshowerE_samedPhi_; - preshowerY_samedPhi_ = preshowerY_samedPhi_/preshowerE_samedPhi_; - preshowerZ_samedPhi_ = preshowerZ_samedPhi_/preshowerE_samedPhi_; - }//if(preshowerNhits_samedPhi_>0) - - - if(preshowerNhits_samedR_>0){ - preshowerX_samedR_ = preshowerX_samedR_/preshowerE_samedR_; - preshowerY_samedR_ = preshowerY_samedR_/preshowerE_samedR_; - preshowerZ_samedR_ = preshowerZ_samedR_/preshowerE_samedR_; - }//if(preshowerNhits_samedR_>0) -} + } ///for(int ih=0; ih 0) { + preshowerX_samedPhi_ = preshowerX_samedPhi_ / preshowerE_samedPhi_; + preshowerY_samedPhi_ = preshowerY_samedPhi_ / preshowerE_samedPhi_; + preshowerZ_samedPhi_ = preshowerZ_samedPhi_ / preshowerE_samedPhi_; + } //if(preshowerNhits_samedPhi_>0) + + if (preshowerNhits_samedR_ > 0) { + preshowerX_samedR_ = preshowerX_samedR_ / preshowerE_samedR_; + preshowerY_samedR_ = preshowerY_samedR_ / preshowerE_samedR_; + preshowerZ_samedR_ = preshowerZ_samedR_ / preshowerE_samedR_; + } //if(preshowerNhits_samedR_>0) +} -double PhotonMVABasedHaloTagger::calAngleBetweenEEAndSubDet(int nhits, double subdetClusX, double subdetClusY, double subdetClusZ){ - +double PhotonMVABasedHaloTagger::calAngleBetweenEEAndSubDet(int nhits, + double subdetClusX, + double subdetClusY, + double subdetClusZ) { ////get the angle of the line joining the ECAL cluster and the subdetector wrt Z axis for any hypothesis - + double angle = -999; - - if(ecalClusNhits_>0 && nhits>0 ){ - - double dR = sqrt( pow(subdetClusX-ecalClusX_,2) + pow(subdetClusY-ecalClusY_,2) + pow(subdetClusZ-ecalClusZ_,2) ); - - double cosTheta = fabs(subdetClusZ-ecalClusZ_)/dR; - + + if (ecalClusNhits_ > 0 && nhits > 0) { + double dR = + sqrt(pow(subdetClusX - ecalClusX_, 2) + pow(subdetClusY - ecalClusY_, 2) + pow(subdetClusZ - ecalClusZ_, 2)); + + double cosTheta = fabs(subdetClusZ - ecalClusZ_) / dR; + angle = acos(cosTheta); } - + return angle; - } - - - From f7d5e90cf92a7becf94e36b145479b8f4bc63931 Mon Sep 17 00:00:00 2001 From: Shilpi Date: Mon, 7 Feb 2022 17:47:52 +0100 Subject: [PATCH 04/13] removed unwanted lines --- .../PhotonIdentification/src/PhotonMVABasedHaloTagger.cc | 9 +-------- 1 file changed, 1 insertion(+), 8 deletions(-) diff --git a/RecoEgamma/PhotonIdentification/src/PhotonMVABasedHaloTagger.cc b/RecoEgamma/PhotonIdentification/src/PhotonMVABasedHaloTagger.cc index 29fdc57e3c23f..f04f8144a9739 100644 --- a/RecoEgamma/PhotonIdentification/src/PhotonMVABasedHaloTagger.cc +++ b/RecoEgamma/PhotonIdentification/src/PhotonMVABasedHaloTagger.cc @@ -154,7 +154,7 @@ void PhotonMVABasedHaloTagger::calphoClusCoordinECAL(const CaloGeometry* geo, double rhY = rechitPoint.y(); double rhZ = rechitPoint.z(); - if ((thresholds == nullptr)) { + if (thresholds == nullptr) { throw cms::Exception("EmptyPFRechHitThresCollection") << "In PhotonMVABasedHaloTagger::calphoClusCoordinECAL, EcalPFRecHitThresholds cannot be = nulptr"; } @@ -323,7 +323,6 @@ void PhotonMVABasedHaloTagger::calmatchedESCoordForBothHypothesis( const GlobalPoint& rechitPoint = geo->getPosition(esrechit->id()); double rhEta = rechitPoint.eta(); - //double rhPhi = rechitPoint.phi(); double rhX = rechitPoint.x(); double rhY = rechitPoint.y(); double rhZ = rechitPoint.z(); @@ -335,11 +334,6 @@ void PhotonMVABasedHaloTagger::calmatchedESCoordForBothHypothesis( if (rhE < noiseThrES_) continue; - /*double rho = sqrt(pow(rhX,2) + pow(rhY,2)); - double dEta = fabs(phoSCEta - rhEta); - double dPhi = deltaPhi(phoSCPhi, rhPhi); - */ - ////try to include RH according to the strips, 11 in X and 11 in Y /////////First calculate RH nearest in phi and eta to that of the photon SC @@ -379,7 +373,6 @@ void PhotonMVABasedHaloTagger::calmatchedESCoordForBothHypothesis( const GlobalPoint& rechitPoint = geo->getPosition(esrechit->id()); double rhEta = rechitPoint.eta(); - double rhPhi = rechitPoint.phi(); double rhX = rechitPoint.x(); double rhY = rechitPoint.y(); double rhZ = rechitPoint.z(); From bb85de93f1236527aab8214040672545faba5a66 Mon Sep 17 00:00:00 2001 From: Shilpi Date: Mon, 7 Feb 2022 22:50:57 +0100 Subject: [PATCH 05/13] missing directory structure --- RecoEgamma/PhotonIdentification/python/mvaHaloVariable_cfi.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/RecoEgamma/PhotonIdentification/python/mvaHaloVariable_cfi.py b/RecoEgamma/PhotonIdentification/python/mvaHaloVariable_cfi.py index 49b99cf575a9a..3067eca8281e7 100644 --- a/RecoEgamma/PhotonIdentification/python/mvaHaloVariable_cfi.py +++ b/RecoEgamma/PhotonIdentification/python/mvaHaloVariable_cfi.py @@ -1,7 +1,7 @@ import FWCore.ParameterSet.Config as cms from RecoEgamma.EgammaIsolationAlgos.egammaHBHERecHitThreshold_cff import egammaHBHERecHit -pathToHaloMVATrainingFile = "RecoEgamma/PhotonIdentification/data/xgboostToTMVA_BHtagger.xml" +pathToHaloMVATrainingFile = "RecoEgamma/PhotonIdentification/data/beamHaloTaggerID/xgboostToTMVA_BHtagger.xml" mvaHaloVariable = cms.PSet( #required inputs trainingFileName = cms.string(pathToHaloMVATrainingFile), From 65d0aa5893425a9ade2fe20c3481a40863c8ef07 Mon Sep 17 00:00:00 2001 From: Shilpi Date: Thu, 17 Feb 2022 21:22:58 +0100 Subject: [PATCH 06/13] Comments from the PR --- .../python/photons_cfi.py | 2 +- .../src/GEDPhotonProducer.cc | 13 +- .../src/PhotonProducer.cc | 11 -- .../interface/PhotonMVABasedHaloTagger.h | 17 +-- .../python/mvaHaloVariable_cfi.py | 2 +- .../src/PhotonMVABasedHaloTagger.cc | 130 ++++++++---------- 6 files changed, 78 insertions(+), 97 deletions(-) diff --git a/RecoEgamma/EgammaPhotonProducers/python/photons_cfi.py b/RecoEgamma/EgammaPhotonProducers/python/photons_cfi.py index 99411cf896367..615aa535e9607 100644 --- a/RecoEgamma/EgammaPhotonProducers/python/photons_cfi.py +++ b/RecoEgamma/EgammaPhotonProducers/python/photons_cfi.py @@ -44,7 +44,7 @@ endcapEcalHits = cms.InputTag("ecalRecHit","EcalRecHitsEE"), preshowerHits = cms.InputTag("ecalPreshowerRecHit","EcalRecHitsES"), runMIPTagger = cms.bool(True), - runMVABasedHaloTagger = cms.bool(True), + runMVABasedHaloTagger = cms.bool(False), highEt = cms.double(100.), minR9Barrel = cms.double(0.94), minR9Endcap = cms.double(0.95), diff --git a/RecoEgamma/EgammaPhotonProducers/src/GEDPhotonProducer.cc b/RecoEgamma/EgammaPhotonProducers/src/GEDPhotonProducer.cc index b85643ac35874..a22e4ed5e814d 100644 --- a/RecoEgamma/EgammaPhotonProducers/src/GEDPhotonProducer.cc +++ b/RecoEgamma/EgammaPhotonProducers/src/GEDPhotonProducer.cc @@ -51,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,9 +69,18 @@ class CacheData { config.outputDim = pset_dnn.getParameter("outputDim"); const auto useEBModelInGap = pset_dnn.getParameter("useEBModelInGap"); photonDNNEstimator = std::make_unique(config, useEBModelInGap); + ///for MVA based beam halo tagger in the EE + } + const auto runMVABasedHaloTagger = conf.getParameter("runMVABasedHaloTagger"); + edm::ParameterSet mvaBasedHaloVariableSet = conf.getParameter("mvaBasedHaloVariableSet"); + auto trainingFileName_ = mvaBasedHaloVariableSet.getParameter("trainingFileName"); + if (runMVABasedHaloTagger) { + TFile* regressionFile = TFile::Open(trainingFileName_.c_str()); + haloTaggerGBR = std::unique_ptr(regressionFile->Get("gbrForest")); } } std::unique_ptr photonDNNEstimator; + std::unique_ptr haloTaggerGBR; }; class GEDPhotonProducer : public edm::stream::EDProducer> { @@ -1069,7 +1080,7 @@ void GEDPhotonProducer::fillPhotonCollection(edm::Event& evt, iSC++; if (runMVABasedHaloTagger_) { ///sets values only for EE, for EB it always returns 1 - float BHmva = photonMVABasedHaloTagger_->calculateMVA(&newCandidate, evt, es); + float BHmva = photonMVABasedHaloTagger_->calculateMVA(&newCandidate, globalCache()->haloTaggerGBR.get(), evt, es); newCandidate.setHaloTaggerMVAVal(BHmva); } diff --git a/RecoEgamma/EgammaPhotonProducers/src/PhotonProducer.cc b/RecoEgamma/EgammaPhotonProducers/src/PhotonProducer.cc index e733bc925bad5..8544ca0f260ca 100644 --- a/RecoEgamma/EgammaPhotonProducers/src/PhotonProducer.cc +++ b/RecoEgamma/EgammaPhotonProducers/src/PhotonProducer.cc @@ -236,12 +236,6 @@ PhotonProducer::PhotonProducer(const edm::ParameterSet& config) edm::ParameterSet mipVariableSet = config.getParameter("mipVariableSet"); photonMIPHaloTagger_.setup(mipVariableSet, consumesCollector()); - if (runMVABasedHaloTagger_) { - edm::ParameterSet mvaBasedHaloVariableSet = config.getParameter("mvaBasedHaloVariableSet"); - photonMVABasedHaloTagger_ = - std::make_unique(mvaBasedHaloVariableSet, consumesCollector()); - } - // Register the product produces(PhotonCollection_); } @@ -498,11 +492,6 @@ void PhotonProducer::fillPhotonCollection(edm::Event& evt, newCandidate.setMIPVariables(mipVar); } - if (runMVABasedHaloTagger_) { ///sets values only for EE, for EB it always returns 1 - double BHmva = photonMVABasedHaloTagger_->calculateMVA(&newCandidate, evt, es); - newCandidate.setHaloTaggerMVAVal(BHmva); - } - /// Pre-selection loose isolation cuts bool isLooseEM = true; if (newCandidate.pt() < highEt_) { diff --git a/RecoEgamma/PhotonIdentification/interface/PhotonMVABasedHaloTagger.h b/RecoEgamma/PhotonIdentification/interface/PhotonMVABasedHaloTagger.h index b25d0e2f3e35f..36fe21bb6027c 100644 --- a/RecoEgamma/PhotonIdentification/interface/PhotonMVABasedHaloTagger.h +++ b/RecoEgamma/PhotonIdentification/interface/PhotonMVABasedHaloTagger.h @@ -5,8 +5,8 @@ 2. JetMET POG: https://indico.cern.ch/event/1027614/contributions/4314949/attachments/2224472/3767396/beamHalo_12April.pdf */ -#ifndef PhotonMVABasedHaloTagger_H -#define PhotonMVABasedHaloTagger_H +#ifndef RecoEgamma_PhotonIdentification_PhotonMVABasedHaloTagger_h +#define RecoEgamma_PhotonIdentification_PhotonMVABasedHaloTagger_h #include "FWCore/Framework/interface/Event.h" #include "FWCore/Framework/interface/EventSetup.h" @@ -27,22 +27,24 @@ class PhotonMVABasedHaloTagger { public: PhotonMVABasedHaloTagger(const edm::ParameterSet& conf, edm::ConsumesCollector&& iC); - virtual ~PhotonMVABasedHaloTagger(){}; - double calculateMVA(const reco::Photon* pho, const edm::Event& iEvent, const edm::EventSetup& es); + double calculateMVA(const reco::Photon* pho, + const GBRForest* gbr_, + const edm::Event& iEvent, + const edm::EventSetup& es); void calphoClusCoordinECAL(const CaloGeometry* geo, const reco::Photon*, const EcalPFRecHitThresholds* thresholds, - edm::Handle ecalrechitCollHandle); + const EcalRecHitCollection& ecalRecHits); void calmatchedHBHECoordForBothHypothesis(const CaloGeometry* geo, const reco::Photon*, - edm::Handle hbheHitHandle); + const HBHERecHitCollection& HBHERecHits); void calmatchedESCoordForBothHypothesis(const CaloGeometry* geo, const reco::Photon*, - edm::Handle esrechitCollHandle); + const EcalRecHitCollection& ESRecHits); double calAngleBetweenEEAndSubDet(int nhits, double subdetClusX, double subdetClusY, double subdetClusZ); @@ -74,7 +76,6 @@ class PhotonMVABasedHaloTagger { edm::EDGetTokenT EEecalCollection_; edm::EDGetTokenT ESCollection_; edm::EDGetTokenT HBHERecHitsCollection_; - std::unique_ptr gbr_; }; #endif // PhotonMVABasedHaloTagger_H diff --git a/RecoEgamma/PhotonIdentification/python/mvaHaloVariable_cfi.py b/RecoEgamma/PhotonIdentification/python/mvaHaloVariable_cfi.py index 3067eca8281e7..a35cd65800713 100644 --- a/RecoEgamma/PhotonIdentification/python/mvaHaloVariable_cfi.py +++ b/RecoEgamma/PhotonIdentification/python/mvaHaloVariable_cfi.py @@ -1,7 +1,7 @@ import FWCore.ParameterSet.Config as cms from RecoEgamma.EgammaIsolationAlgos.egammaHBHERecHitThreshold_cff import egammaHBHERecHit -pathToHaloMVATrainingFile = "RecoEgamma/PhotonIdentification/data/beamHaloTaggerID/xgboostToTMVA_BHtagger.xml" +pathToHaloMVATrainingFile = "$CMSSW_BASE/src/RecoEgamma/PhotonIdentification/data/beamHaloTaggerID/xgboostToTMVA_BHtagger.root" mvaHaloVariable = cms.PSet( #required inputs trainingFileName = cms.string(pathToHaloMVATrainingFile), diff --git a/RecoEgamma/PhotonIdentification/src/PhotonMVABasedHaloTagger.cc b/RecoEgamma/PhotonIdentification/src/PhotonMVABasedHaloTagger.cc index f04f8144a9739..773f2e0912624 100644 --- a/RecoEgamma/PhotonIdentification/src/PhotonMVABasedHaloTagger.cc +++ b/RecoEgamma/PhotonIdentification/src/PhotonMVABasedHaloTagger.cc @@ -30,10 +30,10 @@ PhotonMVABasedHaloTagger::PhotonMVABasedHaloTagger(const edm::ParameterSet& conf recHitEThresholdHE_ = conf.getParameter("recHitEThresholdHE"); noiseThrES_ = conf.getParameter("noiseThrES"); - gbr_ = createGBRForest(trainingFileName_); } double PhotonMVABasedHaloTagger::calculateMVA(const reco::Photon* pho, + const GBRForest* gbr_, const edm::Event& iEvent, const edm::EventSetup& es) { bool isEB = pho->isEB(); @@ -42,23 +42,13 @@ double PhotonMVABasedHaloTagger::calculateMVA(const reco::Photon* pho, return 1.0; /// this MVA is useful and trained only for the EE photons. For EB, there are a lot of other useful handles which can reject beam halo efficiently //rho handle - edm::Handle rhoHandle; - iEvent.getByToken(rhoLabel_, rhoHandle); - - double rho_ = *(rhoHandle.product()); + double rho_ = iEvent.get(rhoLabel_); // Get all the RecHits - edm::Handle barrelHitHandle; - iEvent.getByToken(EBecalCollection_, barrelHitHandle); - - edm::Handle endcapHitHandle; - iEvent.getByToken(EEecalCollection_, endcapHitHandle); - - edm::Handle esHitHandle; - iEvent.getByToken(ESCollection_, esHitHandle); - - edm::Handle hbheHitHandle; - iEvent.getByToken(HBHERecHitsCollection_, hbheHitHandle); + const auto& ecalRecHitsBarrel = iEvent.get(EBecalCollection_); + const auto& ecalRecHitsEndcap = iEvent.get(EEecalCollection_); + const auto& esRecHits = iEvent.get(ESCollection_); + const auto& hbheRecHits = iEvent.get(HBHERecHitsCollection_); //gets geometry pG_ = es.getHandle(geometryToken_); @@ -72,13 +62,13 @@ double PhotonMVABasedHaloTagger::calculateMVA(const reco::Photon* pho, ///calculate the energy weighted X, Y and Z position of the photon cluster if (isEB) - calphoClusCoordinECAL(geo, pho, &thresholds, barrelHitHandle); + calphoClusCoordinECAL(geo, pho, &thresholds, ecalRecHitsBarrel); else - calphoClusCoordinECAL(geo, pho, &thresholds, endcapHitHandle); + calphoClusCoordinECAL(geo, pho, &thresholds, ecalRecHitsEndcap); ///calculate the HBHE cluster position hypothesis - calmatchedHBHECoordForBothHypothesis(geo, pho, hbheHitHandle); - calmatchedESCoordForBothHypothesis(geo, pho, esHitHandle); + calmatchedHBHECoordForBothHypothesis(geo, pho, hbheRecHits); + calmatchedESCoordForBothHypothesis(geo, pho, esRecHits); ///this function works for EE only. Above ones work for EB as well in case later one wants to put a similar function for EB without returning 1 @@ -123,14 +113,14 @@ double PhotonMVABasedHaloTagger::calculateMVA(const reco::Photon* pho, pho->superCluster()->rawEnergy(); vars[14] = rho_; - double BHmva = gbr_->GetAdaBoostClassifier(vars); + double BHmva = gbr_->GetGradBoostClassifier(vars); return BHmva; } void PhotonMVABasedHaloTagger::calphoClusCoordinECAL(const CaloGeometry* geo, const reco::Photon* pho, const EcalPFRecHitThresholds* thresholds, - edm::Handle ecalrechitCollHandle) { + const EcalRecHitCollection& ecalRecHits) { ecalClusX_ = 0; ecalClusY_ = 0; ecalClusZ_ = 0; @@ -140,12 +130,9 @@ void PhotonMVABasedHaloTagger::calphoClusCoordinECAL(const CaloGeometry* geo, double phoSCEta = pho->superCluster()->eta(); double phoSCPhi = pho->superCluster()->phi(); - const EcalRecHitCollection* ecalRecHits = ecalrechitCollHandle.product(); - - for (EcalRecHitCollection::const_iterator ecalrechit = ecalRecHits->begin(); ecalrechit != ecalRecHits->end(); - ecalrechit++) { - auto const det = ecalrechit->id(); - double rhE = ecalrechit->energy(); + for (const auto& ecalrechit : ecalRecHits) { + auto const det = ecalrechit.id(); + double rhE = ecalrechit.energy(); const GlobalPoint& rechitPoint = geo->getPosition(det); double rhEta = rechitPoint.eta(); @@ -160,9 +147,7 @@ void PhotonMVABasedHaloTagger::calphoClusCoordinECAL(const CaloGeometry* geo, } float rhThres = 0.0; - if (thresholds != nullptr) { - rhThres = (*thresholds)[det]; - } + rhThres = (*thresholds)[det]; if (rhE <= rhThres) continue; @@ -170,13 +155,9 @@ void PhotonMVABasedHaloTagger::calphoClusCoordinECAL(const CaloGeometry* geo, if (phoSCEta * rhEta < 0) continue; - //double rho = sqrt(pow(rhX,2) + pow(rhY,2)); - double dPhi = deltaPhi(rhPhi, phoSCPhi); - double dEta = fabs(rhEta - phoSCEta); - - double dR = sqrt(pow(dPhi, 2) + pow(dEta, 2)); + double dR2 = reco::deltaR2(rhEta, rhPhi, phoSCEta, phoSCPhi); - if (dR < 0.2) { + if (dR2 < 0.2 * 0.2) { ecalClusX_ += rhX * rhE; ecalClusY_ += rhY * rhE; ecalClusZ_ += rhZ * rhE; @@ -192,8 +173,9 @@ void PhotonMVABasedHaloTagger::calphoClusCoordinECAL(const CaloGeometry* geo, } //if(ecalClusNhits_>0) } -void PhotonMVABasedHaloTagger::calmatchedHBHECoordForBothHypothesis( - const CaloGeometry* geo, const reco::Photon* pho, edm::Handle hcalhitsCollHBHE) { +void PhotonMVABasedHaloTagger::calmatchedHBHECoordForBothHypothesis(const CaloGeometry* geo, + const reco::Photon* pho, + const HBHERecHitCollection& HBHERecHits) { hcalClusX_samedPhi_ = 0; hcalClusY_samedPhi_ = 0; hcalClusZ_samedPhi_ = 0; @@ -207,15 +189,12 @@ void PhotonMVABasedHaloTagger::calmatchedHBHECoordForBothHypothesis( hcalClusNhits_samedR_ = 0; hcalClusE_samedR_ = 0; - const HBHERecHitCollection* HBHERecHits = hcalhitsCollHBHE.product(); - double phoSCEta = pho->superCluster()->eta(); double phoSCPhi = pho->superCluster()->phi(); // Loop over HBHERecHit's - for (HBHERecHitCollection::const_iterator hbherechit = HBHERecHits->begin(); hbherechit != HBHERecHits->end(); - hbherechit++) { - HcalDetId det = hbherechit->id(); + for (const auto& hbherechit : HBHERecHits) { + HcalDetId det = hbherechit.id(); const GlobalPoint& rechitPoint = geo->getPosition(det); double rhEta = rechitPoint.eta(); @@ -223,7 +202,7 @@ void PhotonMVABasedHaloTagger::calmatchedHBHECoordForBothHypothesis( double rhX = rechitPoint.x(); double rhY = rechitPoint.y(); double rhZ = rechitPoint.z(); - double rhE = hbherechit->energy(); + double rhE = hbherechit.energy(); int depth = det.depth(); @@ -240,16 +219,16 @@ void PhotonMVABasedHaloTagger::calmatchedHBHECoordForBothHypothesis( if (phoSCEta * rhEta < 0) continue; ///Should be on the same side of Z - double rho = sqrt(pow(rhX, 2) + pow(rhY, 2)); - double dEta = fabs(phoSCEta - rhEta); + double rho2 = pow(rhX, 2) + pow(rhY, 2); double dPhi = deltaPhi(phoSCPhi, rhPhi); bool isRHBehindECAL = false; - double dRho = sqrt(pow(rhX - ecalClusX_, 2) + pow(rhY - ecalClusY_, 2)); + double dRho2 = pow(rhX - ecalClusX_, 2) + pow(rhY - ecalClusY_, 2); - if (rho >= 31 && rho <= 172 && dRho <= 26 && - fabs(dPhi) < 0.15) { ///only valid for the EE; this is 26 cm; hit within 3x3 of HCAL centered at the EECAL xtal + if (rho2 >= 31 * 31 && rho2 <= 172 * 172 && dRho2 <= 26 * 26 && + std::abs(dPhi) < + 0.15) { ///only valid for the EE; this is 26 cm; hit within 3x3 of HCAL centered at the EECAL xtal hcalClusX_samedPhi_ += rhX * rhE; hcalClusY_samedPhi_ += rhY * rhE; hcalClusZ_samedPhi_ += rhZ * rhE; @@ -259,9 +238,9 @@ void PhotonMVABasedHaloTagger::calmatchedHBHECoordForBothHypothesis( } //if(rho>=31 && rho<=172) - double dR = sqrt(pow(dEta, 2) + pow(dPhi, 2)); + double dR2 = reco::deltaR2(phoSCEta, phoSCPhi, rhEta, rhPhi); - if (dR < 0.15 && !isRHBehindECAL) { ///dont use hits which are just behind the ECAL in the same phi region + if (dR2 < 0.15 * 0.15 && !isRHBehindECAL) { ///dont use hits which are just behind the ECAL in the same phi region hcalClusX_samedR_ += rhX * rhE; hcalClusY_samedR_ += rhY * rhE; hcalClusZ_samedR_ += rhZ * rhE; @@ -284,8 +263,9 @@ void PhotonMVABasedHaloTagger::calmatchedHBHECoordForBothHypothesis( } //if(hcalClusNhits_samedR_>0) } -void PhotonMVABasedHaloTagger::calmatchedESCoordForBothHypothesis( - const CaloGeometry* geo, const reco::Photon* pho, edm::Handle esrechitCollHandle) { +void PhotonMVABasedHaloTagger::calmatchedESCoordForBothHypothesis(const CaloGeometry* geo, + const reco::Photon* pho, + const EcalRecHitCollection& ESRecHits) { preshowerX_samedPhi_ = 0; preshowerY_samedPhi_ = 0; preshowerZ_samedPhi_ = 0; @@ -298,8 +278,6 @@ void PhotonMVABasedHaloTagger::calmatchedESCoordForBothHypothesis( preshowerNhits_samedR_ = 0; preshowerE_samedR_ = 0; - const EcalRecHitCollection* ESRecHits = esrechitCollHandle.product(); - double phoSCEta = pho->superCluster()->eta(); double phoSCPhi = pho->superCluster()->phi(); @@ -319,14 +297,14 @@ void PhotonMVABasedHaloTagger::calmatchedESCoordForBothHypothesis( double cos_phi = cos(phoSCPhi); double sin_phi = sin(phoSCPhi); - for (EcalRecHitCollection::const_iterator esrechit = ESRecHits->begin(); esrechit != ESRecHits->end(); esrechit++) { - const GlobalPoint& rechitPoint = geo->getPosition(esrechit->id()); + for (const auto& esrechit : ESRecHits) { + const GlobalPoint& rechitPoint = geo->getPosition(esrechit.id()); double rhEta = rechitPoint.eta(); double rhX = rechitPoint.x(); double rhY = rechitPoint.y(); double rhZ = rechitPoint.z(); - double rhE = esrechit->energy(); + double rhE = esrechit.energy(); if (phoSCEta * rhEta < 0) continue; @@ -338,12 +316,14 @@ void PhotonMVABasedHaloTagger::calmatchedESCoordForBothHypothesis( /////////First calculate RH nearest in phi and eta to that of the photon SC //////same phi ----> the X and Y should be similar - double dRho = sqrt(pow(rhX - ecalClusX_, 2) + pow(rhY - ecalClusY_, 2)); + double dRho2 = pow(rhX - ecalClusX_, 2) + pow(rhY - ecalClusY_, 2); - if (dRho < tmpDiffdRho && - dRho < 2.2) { ////i.e. hit is required to be within the ----> seems better match with the data compared to 2.47 + if (dRho2 < tmpDiffdRho && + dRho2 < + 2.2 * + 2.2) { ////i.e. hit is required to be within the ----> seems better match with the data compared to 2.47 - tmpDiffdRho = dRho; + tmpDiffdRho = dRho2; matchX_samephi = rhX; matchY_samephi = rhY; foundESRH_samephi = true; @@ -355,28 +335,28 @@ void PhotonMVABasedHaloTagger::calmatchedESCoordForBothHypothesis( double exp_ESX = cos_phi * exp_ESRho; double exp_ESY = sin_phi * exp_ESRho; - double dRho_samedR = sqrt(pow(rhX - exp_ESX, 2) + pow(rhY - exp_ESY, 2)); + double dRho_samedR2 = pow(rhX - exp_ESX, 2) + pow(rhY - exp_ESY, 2); - if (dRho_samedR < tmpDiffdRho_samedR) { - tmpDiffdRho_samedR = dRho_samedR; + if (dRho_samedR2 < tmpDiffdRho_samedR) { + tmpDiffdRho_samedR = dRho_samedR2; matchX_samedR = rhX; matchY_samedR = rhY; foundESRH_samedR = true; } - } /// for( EcalRecHitCollection::const_iterator esrechit = ESRecHits->begin(); esrechit != ESRecHits->end(); esrechit++ ){ + } /// for (const auto& esrechit : ESRecHits) ////Now calculate the sum in +/- 5 strips in X and y around the matched RH //+/5 strips mean = 5*~2mm = +/-10 mm = 1 cm - for (EcalRecHitCollection::const_iterator esrechit = ESRecHits->begin(); esrechit != ESRecHits->end(); esrechit++) { - const GlobalPoint& rechitPoint = geo->getPosition(esrechit->id()); + for (const auto& esrechit : ESRecHits) { + const GlobalPoint& rechitPoint = geo->getPosition(esrechit.id()); double rhEta = rechitPoint.eta(); double rhX = rechitPoint.x(); double rhY = rechitPoint.y(); double rhZ = rechitPoint.z(); - double rhE = esrechit->energy(); + double rhE = esrechit.energy(); if (phoSCEta * rhEta < 0) continue; @@ -386,8 +366,8 @@ void PhotonMVABasedHaloTagger::calmatchedESCoordForBothHypothesis( bool isRHBehindECAL = false; ///same phi - double dX_samephi = fabs(matchX_samephi - rhX); - double dY_samephi = fabs(matchY_samephi - rhY); + double dX_samephi = std::abs(matchX_samephi - rhX); + double dY_samephi = std::abs(matchY_samephi - rhY); if ((dX_samephi < 1 && dY_samephi < 1) && foundESRH_samephi) { preshowerX_samedPhi_ += rhX * rhE; @@ -399,8 +379,8 @@ void PhotonMVABasedHaloTagger::calmatchedESCoordForBothHypothesis( } ///same dR - double dX_samedR = fabs(matchX_samedR - rhX); - double dY_samedR = fabs(matchY_samedR - rhY); + double dX_samedR = std::abs(matchX_samedR - rhX); + double dY_samedR = std::abs(matchY_samedR - rhY); if (!isRHBehindECAL && foundESRH_samedR && (dX_samedR < 1 && dY_samedR < 1)) { preshowerX_samedR_ += rhX * rhE; @@ -436,7 +416,7 @@ double PhotonMVABasedHaloTagger::calAngleBetweenEEAndSubDet(int nhits, double dR = sqrt(pow(subdetClusX - ecalClusX_, 2) + pow(subdetClusY - ecalClusY_, 2) + pow(subdetClusZ - ecalClusZ_, 2)); - double cosTheta = fabs(subdetClusZ - ecalClusZ_) / dR; + double cosTheta = std::abs(subdetClusZ - ecalClusZ_) / dR; angle = acos(cosTheta); } From c62a11d837fdac359ec4f71ccfdf839708f1609c Mon Sep 17 00:00:00 2001 From: Shilpi Date: Fri, 18 Feb 2022 16:41:50 +0100 Subject: [PATCH 07/13] Corrected full file path --- .../src/GEDPhotonProducer.cc | 8 +++- .../python/mvaHaloVariable_cfi.py | 4 +- .../src/PhotonMVABasedHaloTagger.cc | 41 +++++++++---------- 3 files changed, 27 insertions(+), 26 deletions(-) diff --git a/RecoEgamma/EgammaPhotonProducers/src/GEDPhotonProducer.cc b/RecoEgamma/EgammaPhotonProducers/src/GEDPhotonProducer.cc index a22e4ed5e814d..f08a1d514ad20 100644 --- a/RecoEgamma/EgammaPhotonProducers/src/GEDPhotonProducer.cc +++ b/RecoEgamma/EgammaPhotonProducers/src/GEDPhotonProducer.cc @@ -73,14 +73,18 @@ class CacheData { } const auto runMVABasedHaloTagger = conf.getParameter("runMVABasedHaloTagger"); edm::ParameterSet mvaBasedHaloVariableSet = conf.getParameter("mvaBasedHaloVariableSet"); - auto trainingFileName_ = mvaBasedHaloVariableSet.getParameter("trainingFileName"); + auto trainingFileName_ = mvaBasedHaloVariableSet.getParameter(("trainingFileName")).fullPath().c_str(); if (runMVABasedHaloTagger) { - TFile* regressionFile = TFile::Open(trainingFileName_.c_str()); + TFile* regressionFile = TFile::Open(trainingFileName_); + //std::unique_ptr up(gbrForestFile.Get("gbrForest")); + haloTaggerGBR = std::unique_ptr(regressionFile->Get("gbrForest")); } + } std::unique_ptr photonDNNEstimator; std::unique_ptr haloTaggerGBR; + }; class GEDPhotonProducer : public edm::stream::EDProducer> { diff --git a/RecoEgamma/PhotonIdentification/python/mvaHaloVariable_cfi.py b/RecoEgamma/PhotonIdentification/python/mvaHaloVariable_cfi.py index a35cd65800713..b0e1cb278a953 100644 --- a/RecoEgamma/PhotonIdentification/python/mvaHaloVariable_cfi.py +++ b/RecoEgamma/PhotonIdentification/python/mvaHaloVariable_cfi.py @@ -1,10 +1,10 @@ import FWCore.ParameterSet.Config as cms from RecoEgamma.EgammaIsolationAlgos.egammaHBHERecHitThreshold_cff import egammaHBHERecHit -pathToHaloMVATrainingFile = "$CMSSW_BASE/src/RecoEgamma/PhotonIdentification/data/beamHaloTaggerID/xgboostToTMVA_BHtagger.root" +pathToHaloMVATrainingFile = "RecoEgamma/PhotonIdentification/data/beamHaloTaggerID/xgboostToTMVA_BHtagger.root" mvaHaloVariable = cms.PSet( #required inputs - trainingFileName = cms.string(pathToHaloMVATrainingFile), + trainingFileName = cms.FileInPath(pathToHaloMVATrainingFile), rhoLabel = cms.InputTag("fixedGridRhoFastjetAllTmp"), barrelEcalRecHitCollection = cms.InputTag("ecalRecHit","EcalRecHitsEB"), endcapEcalRecHitCollection = cms.InputTag("ecalRecHit","EcalRecHitsEE"), diff --git a/RecoEgamma/PhotonIdentification/src/PhotonMVABasedHaloTagger.cc b/RecoEgamma/PhotonIdentification/src/PhotonMVABasedHaloTagger.cc index 773f2e0912624..683ba8a6e456c 100644 --- a/RecoEgamma/PhotonIdentification/src/PhotonMVABasedHaloTagger.cc +++ b/RecoEgamma/PhotonIdentification/src/PhotonMVABasedHaloTagger.cc @@ -19,7 +19,6 @@ PhotonMVABasedHaloTagger::PhotonMVABasedHaloTagger(const edm::ParameterSet& conf : geometryToken_(iC.esConsumes()), ecalPFRechitThresholdsToken_(iC.esConsumes()), ecalClusterToolsESGetTokens_(std::move(iC)) { - trainingFileName_ = conf.getParameter("trainingFileName"); rhoLabel_ = iC.consumes(conf.getParameter("rhoLabel")); EBecalCollection_ = iC.consumes(conf.getParameter("barrelEcalRecHitCollection")); @@ -33,7 +32,7 @@ PhotonMVABasedHaloTagger::PhotonMVABasedHaloTagger(const edm::ParameterSet& conf } double PhotonMVABasedHaloTagger::calculateMVA(const reco::Photon* pho, - const GBRForest* gbr_, + const GBRForest* gbr_, const edm::Event& iEvent, const edm::EventSetup& es) { bool isEB = pho->isEB(); @@ -45,8 +44,8 @@ double PhotonMVABasedHaloTagger::calculateMVA(const reco::Photon* pho, double rho_ = iEvent.get(rhoLabel_); // Get all the RecHits - const auto& ecalRecHitsBarrel = iEvent.get(EBecalCollection_); - const auto& ecalRecHitsEndcap = iEvent.get(EEecalCollection_); + const auto& ecalRecHitsBarrel = iEvent.get(EBecalCollection_); + const auto& ecalRecHitsEndcap = iEvent.get(EEecalCollection_); const auto& esRecHits = iEvent.get(ESCollection_); const auto& hbheRecHits = iEvent.get(HBHERecHitsCollection_); @@ -60,6 +59,7 @@ double PhotonMVABasedHaloTagger::calculateMVA(const reco::Photon* pho, noZS::EcalClusterLazyTools lazyToolnoZS( iEvent, ecalClusterToolsESGetTokens_.get(es), EBecalCollection_, EEecalCollection_); + ///calculate the energy weighted X, Y and Z position of the photon cluster if (isEB) calphoClusCoordinECAL(geo, pho, &thresholds, ecalRecHitsBarrel); @@ -130,7 +130,7 @@ void PhotonMVABasedHaloTagger::calphoClusCoordinECAL(const CaloGeometry* geo, double phoSCEta = pho->superCluster()->eta(); double phoSCPhi = pho->superCluster()->phi(); - for (const auto& ecalrechit : ecalRecHits) { + for (const auto& ecalrechit : ecalRecHits){ auto const det = ecalrechit.id(); double rhE = ecalrechit.energy(); const GlobalPoint& rechitPoint = geo->getPosition(det); @@ -149,6 +149,7 @@ void PhotonMVABasedHaloTagger::calphoClusCoordinECAL(const CaloGeometry* geo, float rhThres = 0.0; rhThres = (*thresholds)[det]; + if (rhE <= rhThres) continue; @@ -156,8 +157,9 @@ void PhotonMVABasedHaloTagger::calphoClusCoordinECAL(const CaloGeometry* geo, continue; double dR2 = reco::deltaR2(rhEta, rhPhi, phoSCEta, phoSCPhi); + - if (dR2 < 0.2 * 0.2) { + if (dR2 < 0.2*0.2) { ecalClusX_ += rhX * rhE; ecalClusY_ += rhY * rhE; ecalClusZ_ += rhZ * rhE; @@ -173,9 +175,8 @@ void PhotonMVABasedHaloTagger::calphoClusCoordinECAL(const CaloGeometry* geo, } //if(ecalClusNhits_>0) } -void PhotonMVABasedHaloTagger::calmatchedHBHECoordForBothHypothesis(const CaloGeometry* geo, - const reco::Photon* pho, - const HBHERecHitCollection& HBHERecHits) { +void PhotonMVABasedHaloTagger::calmatchedHBHECoordForBothHypothesis( + const CaloGeometry* geo, const reco::Photon* pho, const HBHERecHitCollection& HBHERecHits) { hcalClusX_samedPhi_ = 0; hcalClusY_samedPhi_ = 0; hcalClusZ_samedPhi_ = 0; @@ -193,7 +194,7 @@ void PhotonMVABasedHaloTagger::calmatchedHBHECoordForBothHypothesis(const CaloGe double phoSCPhi = pho->superCluster()->phi(); // Loop over HBHERecHit's - for (const auto& hbherechit : HBHERecHits) { + for (const auto& hbherechit : HBHERecHits){ HcalDetId det = hbherechit.id(); const GlobalPoint& rechitPoint = geo->getPosition(det); @@ -226,9 +227,8 @@ void PhotonMVABasedHaloTagger::calmatchedHBHECoordForBothHypothesis(const CaloGe double dRho2 = pow(rhX - ecalClusX_, 2) + pow(rhY - ecalClusY_, 2); - if (rho2 >= 31 * 31 && rho2 <= 172 * 172 && dRho2 <= 26 * 26 && - std::abs(dPhi) < - 0.15) { ///only valid for the EE; this is 26 cm; hit within 3x3 of HCAL centered at the EECAL xtal + if (rho2 >= 31*31 && rho2 <= 172*172 && dRho2 <= 26*26 && + std::abs(dPhi) < 0.15) { ///only valid for the EE; this is 26 cm; hit within 3x3 of HCAL centered at the EECAL xtal hcalClusX_samedPhi_ += rhX * rhE; hcalClusY_samedPhi_ += rhY * rhE; hcalClusZ_samedPhi_ += rhZ * rhE; @@ -240,7 +240,7 @@ void PhotonMVABasedHaloTagger::calmatchedHBHECoordForBothHypothesis(const CaloGe double dR2 = reco::deltaR2(phoSCEta, phoSCPhi, rhEta, rhPhi); - if (dR2 < 0.15 * 0.15 && !isRHBehindECAL) { ///dont use hits which are just behind the ECAL in the same phi region + if (dR2 < 0.15*0.15 && !isRHBehindECAL) { ///dont use hits which are just behind the ECAL in the same phi region hcalClusX_samedR_ += rhX * rhE; hcalClusY_samedR_ += rhY * rhE; hcalClusZ_samedR_ += rhZ * rhE; @@ -263,9 +263,8 @@ void PhotonMVABasedHaloTagger::calmatchedHBHECoordForBothHypothesis(const CaloGe } //if(hcalClusNhits_samedR_>0) } -void PhotonMVABasedHaloTagger::calmatchedESCoordForBothHypothesis(const CaloGeometry* geo, - const reco::Photon* pho, - const EcalRecHitCollection& ESRecHits) { +void PhotonMVABasedHaloTagger::calmatchedESCoordForBothHypothesis( + const CaloGeometry* geo, const reco::Photon* pho, const EcalRecHitCollection& ESRecHits) { preshowerX_samedPhi_ = 0; preshowerY_samedPhi_ = 0; preshowerZ_samedPhi_ = 0; @@ -297,7 +296,7 @@ void PhotonMVABasedHaloTagger::calmatchedESCoordForBothHypothesis(const CaloGeom double cos_phi = cos(phoSCPhi); double sin_phi = sin(phoSCPhi); - for (const auto& esrechit : ESRecHits) { + for (const auto& esrechit : ESRecHits){ const GlobalPoint& rechitPoint = geo->getPosition(esrechit.id()); double rhEta = rechitPoint.eta(); @@ -319,9 +318,7 @@ void PhotonMVABasedHaloTagger::calmatchedESCoordForBothHypothesis(const CaloGeom double dRho2 = pow(rhX - ecalClusX_, 2) + pow(rhY - ecalClusY_, 2); if (dRho2 < tmpDiffdRho && - dRho2 < - 2.2 * - 2.2) { ////i.e. hit is required to be within the ----> seems better match with the data compared to 2.47 + dRho2 < 2.2*2.2) { ////i.e. hit is required to be within the ----> seems better match with the data compared to 2.47 tmpDiffdRho = dRho2; matchX_samephi = rhX; @@ -349,7 +346,7 @@ void PhotonMVABasedHaloTagger::calmatchedESCoordForBothHypothesis(const CaloGeom ////Now calculate the sum in +/- 5 strips in X and y around the matched RH //+/5 strips mean = 5*~2mm = +/-10 mm = 1 cm - for (const auto& esrechit : ESRecHits) { + for (const auto& esrechit : ESRecHits){ const GlobalPoint& rechitPoint = geo->getPosition(esrechit.id()); double rhEta = rechitPoint.eta(); From 85d62c58089d3c4aad5b8a5cfffc51893365240e Mon Sep 17 00:00:00 2001 From: Shilpi Date: Fri, 18 Feb 2022 17:25:28 +0100 Subject: [PATCH 08/13] code format checks --- .../src/GEDPhotonProducer.cc | 5 +-- .../src/PhotonMVABasedHaloTagger.cc | 40 ++++++++++--------- 2 files changed, 23 insertions(+), 22 deletions(-) diff --git a/RecoEgamma/EgammaPhotonProducers/src/GEDPhotonProducer.cc b/RecoEgamma/EgammaPhotonProducers/src/GEDPhotonProducer.cc index f08a1d514ad20..490a3029f9147 100644 --- a/RecoEgamma/EgammaPhotonProducers/src/GEDPhotonProducer.cc +++ b/RecoEgamma/EgammaPhotonProducers/src/GEDPhotonProducer.cc @@ -73,18 +73,17 @@ class CacheData { } const auto runMVABasedHaloTagger = conf.getParameter("runMVABasedHaloTagger"); edm::ParameterSet mvaBasedHaloVariableSet = conf.getParameter("mvaBasedHaloVariableSet"); - auto trainingFileName_ = mvaBasedHaloVariableSet.getParameter(("trainingFileName")).fullPath().c_str(); + auto trainingFileName_ = + mvaBasedHaloVariableSet.getParameter(("trainingFileName")).fullPath().c_str(); if (runMVABasedHaloTagger) { TFile* regressionFile = TFile::Open(trainingFileName_); //std::unique_ptr up(gbrForestFile.Get("gbrForest")); haloTaggerGBR = std::unique_ptr(regressionFile->Get("gbrForest")); } - } std::unique_ptr photonDNNEstimator; std::unique_ptr haloTaggerGBR; - }; class GEDPhotonProducer : public edm::stream::EDProducer> { diff --git a/RecoEgamma/PhotonIdentification/src/PhotonMVABasedHaloTagger.cc b/RecoEgamma/PhotonIdentification/src/PhotonMVABasedHaloTagger.cc index 683ba8a6e456c..87b768fe7c199 100644 --- a/RecoEgamma/PhotonIdentification/src/PhotonMVABasedHaloTagger.cc +++ b/RecoEgamma/PhotonIdentification/src/PhotonMVABasedHaloTagger.cc @@ -32,7 +32,7 @@ PhotonMVABasedHaloTagger::PhotonMVABasedHaloTagger(const edm::ParameterSet& conf } double PhotonMVABasedHaloTagger::calculateMVA(const reco::Photon* pho, - const GBRForest* gbr_, + const GBRForest* gbr_, const edm::Event& iEvent, const edm::EventSetup& es) { bool isEB = pho->isEB(); @@ -44,8 +44,8 @@ double PhotonMVABasedHaloTagger::calculateMVA(const reco::Photon* pho, double rho_ = iEvent.get(rhoLabel_); // Get all the RecHits - const auto& ecalRecHitsBarrel = iEvent.get(EBecalCollection_); - const auto& ecalRecHitsEndcap = iEvent.get(EEecalCollection_); + const auto& ecalRecHitsBarrel = iEvent.get(EBecalCollection_); + const auto& ecalRecHitsEndcap = iEvent.get(EEecalCollection_); const auto& esRecHits = iEvent.get(ESCollection_); const auto& hbheRecHits = iEvent.get(HBHERecHitsCollection_); @@ -59,7 +59,6 @@ double PhotonMVABasedHaloTagger::calculateMVA(const reco::Photon* pho, noZS::EcalClusterLazyTools lazyToolnoZS( iEvent, ecalClusterToolsESGetTokens_.get(es), EBecalCollection_, EEecalCollection_); - ///calculate the energy weighted X, Y and Z position of the photon cluster if (isEB) calphoClusCoordinECAL(geo, pho, &thresholds, ecalRecHitsBarrel); @@ -130,7 +129,7 @@ void PhotonMVABasedHaloTagger::calphoClusCoordinECAL(const CaloGeometry* geo, double phoSCEta = pho->superCluster()->eta(); double phoSCPhi = pho->superCluster()->phi(); - for (const auto& ecalrechit : ecalRecHits){ + for (const auto& ecalrechit : ecalRecHits) { auto const det = ecalrechit.id(); double rhE = ecalrechit.energy(); const GlobalPoint& rechitPoint = geo->getPosition(det); @@ -149,7 +148,6 @@ void PhotonMVABasedHaloTagger::calphoClusCoordinECAL(const CaloGeometry* geo, float rhThres = 0.0; rhThres = (*thresholds)[det]; - if (rhE <= rhThres) continue; @@ -157,9 +155,8 @@ void PhotonMVABasedHaloTagger::calphoClusCoordinECAL(const CaloGeometry* geo, continue; double dR2 = reco::deltaR2(rhEta, rhPhi, phoSCEta, phoSCPhi); - - if (dR2 < 0.2*0.2) { + if (dR2 < 0.2 * 0.2) { ecalClusX_ += rhX * rhE; ecalClusY_ += rhY * rhE; ecalClusZ_ += rhZ * rhE; @@ -175,8 +172,9 @@ void PhotonMVABasedHaloTagger::calphoClusCoordinECAL(const CaloGeometry* geo, } //if(ecalClusNhits_>0) } -void PhotonMVABasedHaloTagger::calmatchedHBHECoordForBothHypothesis( - const CaloGeometry* geo, const reco::Photon* pho, const HBHERecHitCollection& HBHERecHits) { +void PhotonMVABasedHaloTagger::calmatchedHBHECoordForBothHypothesis(const CaloGeometry* geo, + const reco::Photon* pho, + const HBHERecHitCollection& HBHERecHits) { hcalClusX_samedPhi_ = 0; hcalClusY_samedPhi_ = 0; hcalClusZ_samedPhi_ = 0; @@ -194,7 +192,7 @@ void PhotonMVABasedHaloTagger::calmatchedHBHECoordForBothHypothesis( double phoSCPhi = pho->superCluster()->phi(); // Loop over HBHERecHit's - for (const auto& hbherechit : HBHERecHits){ + for (const auto& hbherechit : HBHERecHits) { HcalDetId det = hbherechit.id(); const GlobalPoint& rechitPoint = geo->getPosition(det); @@ -227,8 +225,9 @@ void PhotonMVABasedHaloTagger::calmatchedHBHECoordForBothHypothesis( double dRho2 = pow(rhX - ecalClusX_, 2) + pow(rhY - ecalClusY_, 2); - if (rho2 >= 31*31 && rho2 <= 172*172 && dRho2 <= 26*26 && - std::abs(dPhi) < 0.15) { ///only valid for the EE; this is 26 cm; hit within 3x3 of HCAL centered at the EECAL xtal + if (rho2 >= 31 * 31 && rho2 <= 172 * 172 && dRho2 <= 26 * 26 && + std::abs(dPhi) < + 0.15) { ///only valid for the EE; this is 26 cm; hit within 3x3 of HCAL centered at the EECAL xtal hcalClusX_samedPhi_ += rhX * rhE; hcalClusY_samedPhi_ += rhY * rhE; hcalClusZ_samedPhi_ += rhZ * rhE; @@ -240,7 +239,7 @@ void PhotonMVABasedHaloTagger::calmatchedHBHECoordForBothHypothesis( double dR2 = reco::deltaR2(phoSCEta, phoSCPhi, rhEta, rhPhi); - if (dR2 < 0.15*0.15 && !isRHBehindECAL) { ///dont use hits which are just behind the ECAL in the same phi region + if (dR2 < 0.15 * 0.15 && !isRHBehindECAL) { ///dont use hits which are just behind the ECAL in the same phi region hcalClusX_samedR_ += rhX * rhE; hcalClusY_samedR_ += rhY * rhE; hcalClusZ_samedR_ += rhZ * rhE; @@ -263,8 +262,9 @@ void PhotonMVABasedHaloTagger::calmatchedHBHECoordForBothHypothesis( } //if(hcalClusNhits_samedR_>0) } -void PhotonMVABasedHaloTagger::calmatchedESCoordForBothHypothesis( - const CaloGeometry* geo, const reco::Photon* pho, const EcalRecHitCollection& ESRecHits) { +void PhotonMVABasedHaloTagger::calmatchedESCoordForBothHypothesis(const CaloGeometry* geo, + const reco::Photon* pho, + const EcalRecHitCollection& ESRecHits) { preshowerX_samedPhi_ = 0; preshowerY_samedPhi_ = 0; preshowerZ_samedPhi_ = 0; @@ -296,7 +296,7 @@ void PhotonMVABasedHaloTagger::calmatchedESCoordForBothHypothesis( double cos_phi = cos(phoSCPhi); double sin_phi = sin(phoSCPhi); - for (const auto& esrechit : ESRecHits){ + for (const auto& esrechit : ESRecHits) { const GlobalPoint& rechitPoint = geo->getPosition(esrechit.id()); double rhEta = rechitPoint.eta(); @@ -318,7 +318,9 @@ void PhotonMVABasedHaloTagger::calmatchedESCoordForBothHypothesis( double dRho2 = pow(rhX - ecalClusX_, 2) + pow(rhY - ecalClusY_, 2); if (dRho2 < tmpDiffdRho && - dRho2 < 2.2*2.2) { ////i.e. hit is required to be within the ----> seems better match with the data compared to 2.47 + dRho2 < + 2.2 * + 2.2) { ////i.e. hit is required to be within the ----> seems better match with the data compared to 2.47 tmpDiffdRho = dRho2; matchX_samephi = rhX; @@ -346,7 +348,7 @@ void PhotonMVABasedHaloTagger::calmatchedESCoordForBothHypothesis( ////Now calculate the sum in +/- 5 strips in X and y around the matched RH //+/5 strips mean = 5*~2mm = +/-10 mm = 1 cm - for (const auto& esrechit : ESRecHits){ + for (const auto& esrechit : ESRecHits) { const GlobalPoint& rechitPoint = geo->getPosition(esrechit.id()); double rhEta = rechitPoint.eta(); From 0f4515d2721ac7cff3414e4382a05cc4b85274df Mon Sep 17 00:00:00 2001 From: Shilpi Date: Sat, 19 Feb 2022 17:01:29 +0100 Subject: [PATCH 09/13] Taking care of Dangling pointer --- RecoEgamma/EgammaPhotonProducers/src/GEDPhotonProducer.cc | 7 ++----- 1 file changed, 2 insertions(+), 5 deletions(-) diff --git a/RecoEgamma/EgammaPhotonProducers/src/GEDPhotonProducer.cc b/RecoEgamma/EgammaPhotonProducers/src/GEDPhotonProducer.cc index 490a3029f9147..64ce3c5f60ad7 100644 --- a/RecoEgamma/EgammaPhotonProducers/src/GEDPhotonProducer.cc +++ b/RecoEgamma/EgammaPhotonProducers/src/GEDPhotonProducer.cc @@ -73,12 +73,9 @@ class CacheData { } const auto runMVABasedHaloTagger = conf.getParameter("runMVABasedHaloTagger"); edm::ParameterSet mvaBasedHaloVariableSet = conf.getParameter("mvaBasedHaloVariableSet"); - auto trainingFileName_ = - mvaBasedHaloVariableSet.getParameter(("trainingFileName")).fullPath().c_str(); + auto trainingFileName_ = mvaBasedHaloVariableSet.getParameter(("trainingFileName")).fullPath(); if (runMVABasedHaloTagger) { - TFile* regressionFile = TFile::Open(trainingFileName_); - //std::unique_ptr up(gbrForestFile.Get("gbrForest")); - + TFile* regressionFile = TFile::Open(trainingFileName_.c_str()); haloTaggerGBR = std::unique_ptr(regressionFile->Get("gbrForest")); } } From 5322e3352c641b3a5e3baded7a1420ece15945b6 Mon Sep 17 00:00:00 2001 From: Shilpi Date: Mon, 21 Feb 2022 11:02:43 +0100 Subject: [PATCH 10/13] Back to using createGBRForests --- RecoEgamma/EgammaPhotonProducers/src/GEDPhotonProducer.cc | 3 +-- 1 file changed, 1 insertion(+), 2 deletions(-) diff --git a/RecoEgamma/EgammaPhotonProducers/src/GEDPhotonProducer.cc b/RecoEgamma/EgammaPhotonProducers/src/GEDPhotonProducer.cc index 64ce3c5f60ad7..48fd0656704c7 100644 --- a/RecoEgamma/EgammaPhotonProducers/src/GEDPhotonProducer.cc +++ b/RecoEgamma/EgammaPhotonProducers/src/GEDPhotonProducer.cc @@ -75,8 +75,7 @@ class CacheData { edm::ParameterSet mvaBasedHaloVariableSet = conf.getParameter("mvaBasedHaloVariableSet"); auto trainingFileName_ = mvaBasedHaloVariableSet.getParameter(("trainingFileName")).fullPath(); if (runMVABasedHaloTagger) { - TFile* regressionFile = TFile::Open(trainingFileName_.c_str()); - haloTaggerGBR = std::unique_ptr(regressionFile->Get("gbrForest")); + haloTaggerGBR = createGBRForest(trainingFileName_); } } std::unique_ptr photonDNNEstimator; From b901548331d2cef31590a62f498a50d91a00a224 Mon Sep 17 00:00:00 2001 From: Shilpi Date: Sun, 27 Feb 2022 17:40:38 +0100 Subject: [PATCH 11/13] second round of comments --- .../interface/PhotonMVABasedHaloTagger.h | 18 ++++++++++--- .../src/PhotonMVABasedHaloTagger.cc | 27 +++++++++---------- 2 files changed, 28 insertions(+), 17 deletions(-) diff --git a/RecoEgamma/PhotonIdentification/interface/PhotonMVABasedHaloTagger.h b/RecoEgamma/PhotonIdentification/interface/PhotonMVABasedHaloTagger.h index 36fe21bb6027c..1942a6d1a6e8a 100644 --- a/RecoEgamma/PhotonIdentification/interface/PhotonMVABasedHaloTagger.h +++ b/RecoEgamma/PhotonIdentification/interface/PhotonMVABasedHaloTagger.h @@ -33,6 +33,8 @@ class PhotonMVABasedHaloTagger { const edm::Event& iEvent, const edm::EventSetup& es); + +private: void calphoClusCoordinECAL(const CaloGeometry* geo, const reco::Photon*, const EcalPFRecHitThresholds* thresholds, @@ -48,7 +50,7 @@ class PhotonMVABasedHaloTagger { double calAngleBetweenEEAndSubDet(int nhits, double subdetClusX, double subdetClusY, double subdetClusZ); -private: + int hcalClusNhits_samedPhi_, hcalClusNhits_samedR_; int ecalClusNhits_, preshowerNhits_samedPhi_, preshowerNhits_samedR_; double hcalClusX_samedPhi_, hcalClusY_samedPhi_, hcalClusZ_samedPhi_, hcalClusX_samedR_, hcalClusY_samedR_, @@ -64,8 +66,6 @@ class PhotonMVABasedHaloTagger { EgammaHcalIsolation::arrayHB recHitEThresholdHB_; EgammaHcalIsolation::arrayHE recHitEThresholdHE_; - std::string trainingFileName_; - const edm::ESGetToken geometryToken_; const edm::ESGetToken ecalPFRechitThresholdsToken_; const EcalClusterLazyTools::ESGetTokens ecalClusterToolsESGetTokens_; @@ -76,6 +76,18 @@ class PhotonMVABasedHaloTagger { edm::EDGetTokenT EEecalCollection_; edm::EDGetTokenT ESCollection_; edm::EDGetTokenT 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 diff --git a/RecoEgamma/PhotonIdentification/src/PhotonMVABasedHaloTagger.cc b/RecoEgamma/PhotonIdentification/src/PhotonMVABasedHaloTagger.cc index 87b768fe7c199..a5a642555aff5 100644 --- a/RecoEgamma/PhotonIdentification/src/PhotonMVABasedHaloTagger.cc +++ b/RecoEgamma/PhotonIdentification/src/PhotonMVABasedHaloTagger.cc @@ -156,7 +156,7 @@ void PhotonMVABasedHaloTagger::calphoClusCoordinECAL(const CaloGeometry* geo, double dR2 = reco::deltaR2(rhEta, rhPhi, phoSCEta, phoSCPhi); - if (dR2 < 0.2 * 0.2) { + if (dR2 < dr2Max_ECALClus_) { ecalClusX_ += rhX * rhE; ecalClusY_ += rhY * rhE; ecalClusZ_ += rhZ * rhE; @@ -206,9 +206,11 @@ void PhotonMVABasedHaloTagger::calmatchedHBHECoordForBothHypothesis(const CaloGe int depth = det.depth(); if ((det.subdet() == HcalBarrel and (depth < 1 or depth > int(recHitEThresholdHB_.size()))) or - (det.subdet() == HcalEndcap and (depth < 1 or depth > int(recHitEThresholdHE_.size())))) + (det.subdet() == HcalEndcap and (depth < 1 or depth > int(recHitEThresholdHE_.size())))){ edm::LogWarning("PhotonMVABasedHaloTagger") - << " hit in subdet " << det.subdet() << " has an unaccounted for depth of " << depth << "!!"; + << " hit in subdet " << det.subdet() << " has an unaccounted for depth of " << depth << "!! Leaving this hit!!"; + continue; + } const bool goodHBe = det.subdet() == HcalBarrel and rhE > recHitEThresholdHB_[depth - 1]; const bool goodHEe = det.subdet() == HcalEndcap and rhE > recHitEThresholdHE_[depth - 1]; @@ -224,10 +226,9 @@ void PhotonMVABasedHaloTagger::calmatchedHBHECoordForBothHypothesis(const CaloGe bool isRHBehindECAL = false; double dRho2 = pow(rhX - ecalClusX_, 2) + pow(rhY - ecalClusY_, 2); - - if (rho2 >= 31 * 31 && rho2 <= 172 * 172 && dRho2 <= 26 * 26 && - std::abs(dPhi) < - 0.15) { ///only valid for the EE; this is 26 cm; hit within 3x3 of HCAL centered at the EECAL xtal + + ///only valid for the EE; this is 26 cm; hit within 3x3 of HCAL centered at the EECAL xtal + if (rho2 >= rho2Min_ECALpos_ && rho2 <= rho2Max_ECALpos_ && dRho2 <= dRho2Max_HCALClus_SamePhi_ && std::abs(dPhi) < dPhiMax_HCALClus_SamePhi_) { hcalClusX_samedPhi_ += rhX * rhE; hcalClusY_samedPhi_ += rhY * rhE; hcalClusZ_samedPhi_ += rhZ * rhE; @@ -239,7 +240,7 @@ void PhotonMVABasedHaloTagger::calmatchedHBHECoordForBothHypothesis(const CaloGe double dR2 = reco::deltaR2(phoSCEta, phoSCPhi, rhEta, rhPhi); - if (dR2 < 0.15 * 0.15 && !isRHBehindECAL) { ///dont use hits which are just behind the ECAL in the same phi region + if (dR2 < dR2Max_HCALClus_SamePhi_ && !isRHBehindECAL) { ///dont use hits which are just behind the ECAL in the same phi region hcalClusX_samedR_ += rhX * rhE; hcalClusY_samedR_ += rhY * rhE; hcalClusZ_samedR_ += rhZ * rhE; @@ -315,12 +316,10 @@ void PhotonMVABasedHaloTagger::calmatchedESCoordForBothHypothesis(const CaloGeom /////////First calculate RH nearest in phi and eta to that of the photon SC //////same phi ----> the X and Y should be similar + ////i.e. hit is required to be within the ----> seems better match with the data compared to 2.47 double dRho2 = pow(rhX - ecalClusX_, 2) + pow(rhY - ecalClusY_, 2); - if (dRho2 < tmpDiffdRho && - dRho2 < - 2.2 * - 2.2) { ////i.e. hit is required to be within the ----> seems better match with the data compared to 2.47 + if (dRho2 < tmpDiffdRho && dRho2 < dRho2Max_ESClus_) { tmpDiffdRho = dRho2; matchX_samephi = rhX; @@ -368,7 +367,7 @@ void PhotonMVABasedHaloTagger::calmatchedESCoordForBothHypothesis(const CaloGeom double dX_samephi = std::abs(matchX_samephi - rhX); double dY_samephi = std::abs(matchY_samephi - rhY); - if ((dX_samephi < 1 && dY_samephi < 1) && foundESRH_samephi) { + if ((dX_samephi < dXY_ESClus_SamePhi_ && dY_samephi < dXY_ESClus_SamePhi_) && foundESRH_samephi) { preshowerX_samedPhi_ += rhX * rhE; preshowerY_samedPhi_ += rhY * rhE; preshowerZ_samedPhi_ += rhZ * rhE; @@ -381,7 +380,7 @@ void PhotonMVABasedHaloTagger::calmatchedESCoordForBothHypothesis(const CaloGeom double dX_samedR = std::abs(matchX_samedR - rhX); double dY_samedR = std::abs(matchY_samedR - rhY); - if (!isRHBehindECAL && foundESRH_samedR && (dX_samedR < 1 && dY_samedR < 1)) { + if (!isRHBehindECAL && foundESRH_samedR && (dX_samedR < dXY_ESClus_SamedR_ && dY_samedR < dXY_ESClus_SamedR_)) { preshowerX_samedR_ += rhX * rhE; preshowerY_samedR_ += rhY * rhE; preshowerZ_samedR_ += rhZ * rhE; From 95736ccffa1305759e6f0b15196fb057d5eff977 Mon Sep 17 00:00:00 2001 From: Shilpi Date: Sun, 27 Feb 2022 17:44:00 +0100 Subject: [PATCH 12/13] second round of comments --- .../src/GEDPhotonProducer.cc | 2 +- .../interface/PhotonMVABasedHaloTagger.h | 17 +++++++---------- .../src/PhotonMVABasedHaloTagger.cc | 14 ++++++++------ 3 files changed, 16 insertions(+), 17 deletions(-) diff --git a/RecoEgamma/EgammaPhotonProducers/src/GEDPhotonProducer.cc b/RecoEgamma/EgammaPhotonProducers/src/GEDPhotonProducer.cc index 48fd0656704c7..e2cad868ccab1 100644 --- a/RecoEgamma/EgammaPhotonProducers/src/GEDPhotonProducer.cc +++ b/RecoEgamma/EgammaPhotonProducers/src/GEDPhotonProducer.cc @@ -69,8 +69,8 @@ class CacheData { config.outputDim = pset_dnn.getParameter("outputDim"); const auto useEBModelInGap = pset_dnn.getParameter("useEBModelInGap"); photonDNNEstimator = std::make_unique(config, useEBModelInGap); - ///for MVA based beam halo tagger in the EE } + ///for MVA based beam halo tagger in the EE const auto runMVABasedHaloTagger = conf.getParameter("runMVABasedHaloTagger"); edm::ParameterSet mvaBasedHaloVariableSet = conf.getParameter("mvaBasedHaloVariableSet"); auto trainingFileName_ = mvaBasedHaloVariableSet.getParameter(("trainingFileName")).fullPath(); diff --git a/RecoEgamma/PhotonIdentification/interface/PhotonMVABasedHaloTagger.h b/RecoEgamma/PhotonIdentification/interface/PhotonMVABasedHaloTagger.h index 1942a6d1a6e8a..935bd9e5fead4 100644 --- a/RecoEgamma/PhotonIdentification/interface/PhotonMVABasedHaloTagger.h +++ b/RecoEgamma/PhotonIdentification/interface/PhotonMVABasedHaloTagger.h @@ -33,7 +33,6 @@ class PhotonMVABasedHaloTagger { const edm::Event& iEvent, const edm::EventSetup& es); - private: void calphoClusCoordinECAL(const CaloGeometry* geo, const reco::Photon*, @@ -50,7 +49,6 @@ class PhotonMVABasedHaloTagger { 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_, @@ -78,16 +76,15 @@ class PhotonMVABasedHaloTagger { edm::EDGetTokenT 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 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 - + 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 diff --git a/RecoEgamma/PhotonIdentification/src/PhotonMVABasedHaloTagger.cc b/RecoEgamma/PhotonIdentification/src/PhotonMVABasedHaloTagger.cc index a5a642555aff5..ce78c50564a19 100644 --- a/RecoEgamma/PhotonIdentification/src/PhotonMVABasedHaloTagger.cc +++ b/RecoEgamma/PhotonIdentification/src/PhotonMVABasedHaloTagger.cc @@ -206,9 +206,10 @@ void PhotonMVABasedHaloTagger::calmatchedHBHECoordForBothHypothesis(const CaloGe int depth = det.depth(); if ((det.subdet() == HcalBarrel and (depth < 1 or depth > int(recHitEThresholdHB_.size()))) or - (det.subdet() == HcalEndcap and (depth < 1 or depth > int(recHitEThresholdHE_.size())))){ + (det.subdet() == HcalEndcap and (depth < 1 or depth > int(recHitEThresholdHE_.size())))) { edm::LogWarning("PhotonMVABasedHaloTagger") - << " hit in subdet " << det.subdet() << " has an unaccounted for depth of " << depth << "!! Leaving this hit!!"; + << " hit in subdet " << det.subdet() << " has an unaccounted for depth of " << depth + << "!! Leaving this hit!!"; continue; } @@ -226,9 +227,10 @@ void PhotonMVABasedHaloTagger::calmatchedHBHECoordForBothHypothesis(const CaloGe bool isRHBehindECAL = false; double dRho2 = pow(rhX - ecalClusX_, 2) + pow(rhY - ecalClusY_, 2); - + ///only valid for the EE; this is 26 cm; hit within 3x3 of HCAL centered at the EECAL xtal - if (rho2 >= rho2Min_ECALpos_ && rho2 <= rho2Max_ECALpos_ && dRho2 <= dRho2Max_HCALClus_SamePhi_ && std::abs(dPhi) < dPhiMax_HCALClus_SamePhi_) { + if (rho2 >= rho2Min_ECALpos_ && rho2 <= rho2Max_ECALpos_ && dRho2 <= dRho2Max_HCALClus_SamePhi_ && + std::abs(dPhi) < dPhiMax_HCALClus_SamePhi_) { hcalClusX_samedPhi_ += rhX * rhE; hcalClusY_samedPhi_ += rhY * rhE; hcalClusZ_samedPhi_ += rhZ * rhE; @@ -240,7 +242,8 @@ void PhotonMVABasedHaloTagger::calmatchedHBHECoordForBothHypothesis(const CaloGe double dR2 = reco::deltaR2(phoSCEta, phoSCPhi, rhEta, rhPhi); - if (dR2 < dR2Max_HCALClus_SamePhi_ && !isRHBehindECAL) { ///dont use hits which are just behind the ECAL in the same phi region + if (dR2 < dR2Max_HCALClus_SamePhi_ && + !isRHBehindECAL) { ///dont use hits which are just behind the ECAL in the same phi region hcalClusX_samedR_ += rhX * rhE; hcalClusY_samedR_ += rhY * rhE; hcalClusZ_samedR_ += rhZ * rhE; @@ -320,7 +323,6 @@ void PhotonMVABasedHaloTagger::calmatchedESCoordForBothHypothesis(const CaloGeom double dRho2 = pow(rhX - ecalClusX_, 2) + pow(rhY - ecalClusY_, 2); if (dRho2 < tmpDiffdRho && dRho2 < dRho2Max_ESClus_) { - tmpDiffdRho = dRho2; matchX_samephi = rhX; matchY_samephi = rhY; From 73220c5a171602299e1b9a676d6afde751651f84 Mon Sep 17 00:00:00 2001 From: Shilpi Date: Wed, 2 Mar 2022 19:13:31 +0100 Subject: [PATCH 13/13] Andrea's comments --- .../src/PhotonMVABasedHaloTagger.cc | 98 ++++++++++--------- 1 file changed, 50 insertions(+), 48 deletions(-) diff --git a/RecoEgamma/PhotonIdentification/src/PhotonMVABasedHaloTagger.cc b/RecoEgamma/PhotonIdentification/src/PhotonMVABasedHaloTagger.cc index ce78c50564a19..afaffacd73b56 100644 --- a/RecoEgamma/PhotonIdentification/src/PhotonMVABasedHaloTagger.cc +++ b/RecoEgamma/PhotonIdentification/src/PhotonMVABasedHaloTagger.cc @@ -145,8 +145,7 @@ void PhotonMVABasedHaloTagger::calphoClusCoordinECAL(const CaloGeometry* geo, << "In PhotonMVABasedHaloTagger::calphoClusCoordinECAL, EcalPFRecHitThresholds cannot be = nulptr"; } - float rhThres = 0.0; - rhThres = (*thresholds)[det]; + float rhThres = (*thresholds)[det]; if (rhE <= rhThres) continue; @@ -221,36 +220,37 @@ void PhotonMVABasedHaloTagger::calmatchedHBHECoordForBothHypothesis(const CaloGe if (phoSCEta * rhEta < 0) continue; ///Should be on the same side of Z - double rho2 = pow(rhX, 2) + pow(rhY, 2); double dPhi = deltaPhi(phoSCPhi, rhPhi); - bool isRHBehindECAL = false; - - double dRho2 = pow(rhX - ecalClusX_, 2) + pow(rhY - ecalClusY_, 2); - ///only valid for the EE; this is 26 cm; hit within 3x3 of HCAL centered at the EECAL xtal - if (rho2 >= rho2Min_ECALpos_ && rho2 <= rho2Max_ECALpos_ && dRho2 <= dRho2Max_HCALClus_SamePhi_ && - std::abs(dPhi) < dPhiMax_HCALClus_SamePhi_) { - hcalClusX_samedPhi_ += rhX * rhE; - hcalClusY_samedPhi_ += rhY * rhE; - hcalClusZ_samedPhi_ += rhZ * rhE; - hcalClusE_samedPhi_ += rhE; - hcalClusNhits_samedPhi_++; - isRHBehindECAL = true; - + bool isRHBehindECAL = std::abs(dPhi) < dPhiMax_HCALClus_SamePhi_; + if (isRHBehindECAL) { + double rho2 = pow(rhX, 2) + pow(rhY, 2); + isRHBehindECAL &= (rho2 >= rho2Min_ECALpos_ && rho2 <= rho2Max_ECALpos_); + if (isRHBehindECAL) { + double dRho2 = pow(rhX - ecalClusX_, 2) + pow(rhY - ecalClusY_, 2); + isRHBehindECAL &= dRho2 <= dRho2Max_HCALClus_SamePhi_; + if (isRHBehindECAL) { + hcalClusX_samedPhi_ += rhX * rhE; + hcalClusY_samedPhi_ += rhY * rhE; + hcalClusZ_samedPhi_ += rhZ * rhE; + hcalClusE_samedPhi_ += rhE; + hcalClusNhits_samedPhi_++; + } + } } //if(rho>=31 && rho<=172) - double dR2 = reco::deltaR2(phoSCEta, phoSCPhi, rhEta, rhPhi); - - if (dR2 < dR2Max_HCALClus_SamePhi_ && - !isRHBehindECAL) { ///dont use hits which are just behind the ECAL in the same phi region - hcalClusX_samedR_ += rhX * rhE; - hcalClusY_samedR_ += rhY * rhE; - hcalClusZ_samedR_ += rhZ * rhE; - hcalClusE_samedR_ += rhE; - hcalClusNhits_samedR_++; + ///dont use hits which are just behind the ECAL in the same phi region + if (!isRHBehindECAL) { + double dR2 = reco::deltaR2(phoSCEta, phoSCPhi, rhEta, rhPhi); + if (dR2 < dR2Max_HCALClus_SamePhi_) { + hcalClusX_samedR_ += rhX * rhE; + hcalClusY_samedR_ += rhY * rhE; + hcalClusZ_samedR_ += rhZ * rhE; + hcalClusE_samedR_ += rhE; + hcalClusNhits_samedR_++; + } } - } //for(int ih=0; ih 0) { @@ -363,34 +363,36 @@ void PhotonMVABasedHaloTagger::calmatchedESCoordForBothHypothesis(const CaloGeom if (rhE < noiseThrES_) continue; - bool isRHBehindECAL = false; - ///same phi - double dX_samephi = std::abs(matchX_samephi - rhX); - double dY_samephi = std::abs(matchY_samephi - rhY); - - if ((dX_samephi < dXY_ESClus_SamePhi_ && dY_samephi < dXY_ESClus_SamePhi_) && foundESRH_samephi) { - preshowerX_samedPhi_ += rhX * rhE; - preshowerY_samedPhi_ += rhY * rhE; - preshowerZ_samedPhi_ += rhZ * rhE; - preshowerE_samedPhi_ += rhE; - preshowerNhits_samedPhi_++; - isRHBehindECAL = true; + bool isRHBehindECAL = foundESRH_samephi; + if (isRHBehindECAL) { + double dX_samephi = std::abs(matchX_samephi - rhX); + double dY_samephi = std::abs(matchY_samephi - rhY); + isRHBehindECAL &= (dX_samephi < dXY_ESClus_SamePhi_ && dY_samephi < dXY_ESClus_SamePhi_); + if (isRHBehindECAL) { + preshowerX_samedPhi_ += rhX * rhE; + preshowerY_samedPhi_ += rhY * rhE; + preshowerZ_samedPhi_ += rhZ * rhE; + preshowerE_samedPhi_ += rhE; + preshowerNhits_samedPhi_++; + } } ///same dR - double dX_samedR = std::abs(matchX_samedR - rhX); - double dY_samedR = std::abs(matchY_samedR - rhY); - - if (!isRHBehindECAL && foundESRH_samedR && (dX_samedR < dXY_ESClus_SamedR_ && dY_samedR < dXY_ESClus_SamedR_)) { - preshowerX_samedR_ += rhX * rhE; - preshowerY_samedR_ += rhY * rhE; - preshowerZ_samedR_ += rhZ * rhE; - preshowerE_samedR_ += rhE; - preshowerNhits_samedR_++; + if (!isRHBehindECAL && foundESRH_samedR) { + double dX_samedR = std::abs(matchX_samedR - rhX); + double dY_samedR = std::abs(matchY_samedR - rhY); + + if (dX_samedR < dXY_ESClus_SamedR_ && dY_samedR < dXY_ESClus_SamedR_) { + preshowerX_samedR_ += rhX * rhE; + preshowerY_samedR_ += rhY * rhE; + preshowerZ_samedR_ += rhZ * rhE; + preshowerE_samedR_ += rhE; + preshowerNhits_samedR_++; + } } - } ///for(int ih=0; ih 0) { preshowerX_samedPhi_ = preshowerX_samedPhi_ / preshowerE_samedPhi_; preshowerY_samedPhi_ = preshowerY_samedPhi_ / preshowerE_samedPhi_;