File indexing completed on 2024-10-08 05:11:59
0001
0002
0003
0004
0005
0006
0007
0008 #include "RecoEgamma/PhotonIdentification/interface/PhotonMVABasedHaloTagger.h"
0009 #include "DataFormats/DetId/interface/DetId.h"
0010 #include "DataFormats/EcalDetId/interface/EBDetId.h"
0011 #include "DataFormats/EcalRecHit/interface/EcalRecHit.h"
0012
0013 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0014 #include "RecoEcal/EgammaCoreTools/interface/EcalClusterTools.h"
0015 #include <DataFormats/Math/interface/deltaPhi.h>
0016 #include "CommonTools/MVAUtils/interface/GBRForestTools.h"
0017
0018 PhotonMVABasedHaloTagger::PhotonMVABasedHaloTagger(const edm::ParameterSet& conf, edm::ConsumesCollector&& iC)
0019 : geometryToken_(iC.esConsumes()), ecalPFRechitThresholdsToken_(iC.esConsumes()), ecalClusterToolsESGetTokens_(iC) {
0020 rhoLabel_ = iC.consumes<double>(conf.getParameter<edm::InputTag>("rhoLabel"));
0021
0022 EBecalCollection_ = iC.consumes<EcalRecHitCollection>(conf.getParameter<edm::InputTag>("barrelEcalRecHitCollection"));
0023 EEecalCollection_ = iC.consumes<EcalRecHitCollection>(conf.getParameter<edm::InputTag>("endcapEcalRecHitCollection"));
0024 ESCollection_ = iC.consumes<EcalRecHitCollection>(conf.getParameter<edm::InputTag>("esRecHitCollection"));
0025 HBHERecHitsCollection_ = iC.consumes<HBHERecHitCollection>(conf.getParameter<edm::InputTag>("HBHERecHitsCollection"));
0026 recHitEThresholdHB_ = conf.getParameter<EgammaHcalIsolation::arrayHB>("recHitEThresholdHB");
0027 recHitEThresholdHE_ = conf.getParameter<EgammaHcalIsolation::arrayHE>("recHitEThresholdHE");
0028
0029 noiseThrES_ = conf.getParameter<double>("noiseThrES");
0030 }
0031
0032 double PhotonMVABasedHaloTagger::calculateMVA(const reco::Photon* pho,
0033 const GBRForest* gbr_,
0034 const edm::Event& iEvent,
0035 const edm::EventSetup& es) const {
0036 bool isEB = pho->isEB();
0037
0038 if (isEB)
0039 return 1.0;
0040
0041
0042 double rho_ = iEvent.get(rhoLabel_);
0043
0044
0045 const auto& ecalRecHitsBarrel = iEvent.get(EBecalCollection_);
0046 const auto& ecalRecHitsEndcap = iEvent.get(EEecalCollection_);
0047 const auto& esRecHits = iEvent.get(ESCollection_);
0048 const auto& hbheRecHits = iEvent.get(HBHERecHitsCollection_);
0049
0050
0051 const CaloGeometry& geo = es.getData(geometryToken_);
0052
0053
0054 auto const& thresholds = es.getData(ecalPFRechitThresholdsToken_);
0055
0056 noZS::EcalClusterLazyTools lazyToolnoZS(
0057 iEvent, ecalClusterToolsESGetTokens_.get(es), EBecalCollection_, EEecalCollection_);
0058
0059
0060 EcalClus ecalClus;
0061 if (isEB)
0062 ecalClus = calphoClusCoordinECAL(geo, pho, &thresholds, ecalRecHitsBarrel);
0063 else
0064 ecalClus = calphoClusCoordinECAL(geo, pho, &thresholds, ecalRecHitsEndcap);
0065
0066
0067 auto hcalClus = calmatchedHBHECoordForBothHypothesis(geo, pho, hbheRecHits, ecalClus);
0068 auto preshower = calmatchedESCoordForBothHypothesis(geo, pho, esRecHits, ecalClus);
0069
0070
0071
0072 double angle_EE_HE_samedPhi = calAngleBetweenEEAndSubDet(
0073 hcalClus.samedPhi_,
0074 ecalClus);
0075
0076 double angle_EE_HE_samedR = calAngleBetweenEEAndSubDet(hcalClus.samedR_, ecalClus);
0077
0078 double angle_EE_ES_samedPhi = calAngleBetweenEEAndSubDet(preshower.samedPhi_, ecalClus);
0079
0080 double angle_EE_ES_samedR = calAngleBetweenEEAndSubDet(preshower.samedR_, ecalClus);
0081
0082
0083
0084 const auto& vCov = lazyToolnoZS.localCovariances(*(pho->superCluster()->seed()));
0085 double spp = (std::isnan(vCov[2]) ? 0. : sqrt(vCov[2]));
0086
0087
0088
0089 float vars[15];
0090
0091 vars[0] = preshower.samedPhi_.e_;
0092 vars[1] = hcalClus.samedPhi_.e_;
0093 vars[2] = preshower.samedR_.e_;
0094 vars[3] = hcalClus.samedR_.e_;
0095 vars[4] = pho->full5x5_r9();
0096 vars[5] = pho->superCluster()->etaWidth();
0097 vars[6] = pho->superCluster()->phiWidth();
0098 vars[7] = pho->full5x5_sigmaIetaIeta();
0099 vars[8] = spp;
0100 vars[9] = angle_EE_ES_samedR;
0101 vars[10] = angle_EE_HE_samedR;
0102 vars[11] = angle_EE_ES_samedPhi;
0103 vars[12] = angle_EE_HE_samedPhi;
0104 vars[13] = (pho->superCluster()->preshowerEnergyPlane1() + pho->superCluster()->preshowerEnergyPlane2()) /
0105 pho->superCluster()->rawEnergy();
0106 vars[14] = rho_;
0107
0108 double BHmva = gbr_->GetGradBoostClassifier(vars);
0109 return BHmva;
0110 }
0111
0112 PhotonMVABasedHaloTagger::EcalClus PhotonMVABasedHaloTagger::calphoClusCoordinECAL(
0113 const CaloGeometry& geo,
0114 const reco::Photon* pho,
0115 const EcalPFRecHitThresholds* thresholds,
0116 const EcalRecHitCollection& ecalRecHits) const {
0117 EcalClus ecalClus;
0118
0119 double phoSCEta = pho->superCluster()->eta();
0120 double phoSCPhi = pho->superCluster()->phi();
0121
0122 if (thresholds == nullptr) {
0123 throw cms::Exception("EmptyPFRechHitThresCollection")
0124 << "In PhotonMVABasedHaloTagger::calphoClusCoordinECAL, EcalPFRecHitThresholds cannot be = nulptr";
0125 }
0126
0127 for (const auto& ecalrechit : ecalRecHits) {
0128 auto const det = ecalrechit.id();
0129 double rhE = ecalrechit.energy();
0130 float rhThres = (*thresholds)[det];
0131 if (rhE <= rhThres)
0132 continue;
0133
0134 const GlobalPoint& rechitPoint = geo.getPosition(det);
0135
0136 double rhEta = rechitPoint.eta();
0137 double rhPhi = rechitPoint.phi();
0138 double rhX = rechitPoint.x();
0139 double rhY = rechitPoint.y();
0140 double rhZ = rechitPoint.z();
0141
0142 if (phoSCEta * rhEta < 0)
0143 continue;
0144
0145 double dR2 = reco::deltaR2(rhEta, rhPhi, phoSCEta, phoSCPhi);
0146
0147 if (dR2 < dr2Max_ECALClus_) {
0148 ecalClus.x_ += rhX * rhE;
0149 ecalClus.y_ += rhY * rhE;
0150 ecalClus.z_ += rhZ * rhE;
0151 ecalClus.e_ += rhE;
0152 ecalClus.nHits_++;
0153 }
0154 }
0155
0156 if (ecalClus.nHits_ > 0) {
0157 ecalClus.x_ = ecalClus.x_ / ecalClus.e_;
0158 ecalClus.y_ = ecalClus.y_ / ecalClus.e_;
0159 ecalClus.z_ = ecalClus.z_ / ecalClus.e_;
0160 }
0161
0162 return ecalClus;
0163 }
0164
0165 PhotonMVABasedHaloTagger::HcalHyp PhotonMVABasedHaloTagger::calmatchedHBHECoordForBothHypothesis(
0166 const CaloGeometry& geo,
0167 const reco::Photon* pho,
0168 const HBHERecHitCollection& HBHERecHits,
0169 const EcalClus& ecalClus) const {
0170 HcalHyp returnValue;
0171
0172 double phoSCEta = pho->superCluster()->eta();
0173 double phoSCPhi = pho->superCluster()->phi();
0174
0175
0176 for (const auto& hbherechit : HBHERecHits) {
0177 HcalDetId det = hbherechit.id();
0178 const GlobalPoint& rechitPoint = geo.getPosition(det);
0179
0180 double rhEta = rechitPoint.eta();
0181 double rhPhi = rechitPoint.phi();
0182 double rhX = rechitPoint.x();
0183 double rhY = rechitPoint.y();
0184 double rhZ = rechitPoint.z();
0185 double rhE = hbherechit.energy();
0186
0187 int depth = det.depth();
0188
0189 if ((det.subdet() == HcalBarrel and (depth < 1 or depth > int(recHitEThresholdHB_.size()))) or
0190 (det.subdet() == HcalEndcap and (depth < 1 or depth > int(recHitEThresholdHE_.size())))) {
0191 edm::LogWarning("PhotonMVABasedHaloTagger")
0192 << " hit in subdet " << det.subdet() << " has an unaccounted for depth of " << depth
0193 << "!! Leaving this hit!!";
0194 continue;
0195 }
0196
0197 const bool goodHBe = det.subdet() == HcalBarrel and rhE > recHitEThresholdHB_[depth - 1];
0198 const bool goodHEe = det.subdet() == HcalEndcap and rhE > recHitEThresholdHE_[depth - 1];
0199 if (!(goodHBe or goodHEe))
0200 continue;
0201
0202 if (phoSCEta * rhEta < 0)
0203 continue;
0204
0205 double dPhi = deltaPhi(phoSCPhi, rhPhi);
0206
0207
0208 bool isRHBehindECAL = std::abs(dPhi) < dPhiMax_HCALClus_SamePhi_;
0209 if (isRHBehindECAL) {
0210 double rho2 = pow(rhX, 2) + pow(rhY, 2);
0211 isRHBehindECAL &= (rho2 >= rho2Min_ECALpos_ && rho2 <= rho2Max_ECALpos_);
0212 if (isRHBehindECAL) {
0213 double dRho2 = pow(rhX - ecalClus.x_, 2) + pow(rhY - ecalClus.y_, 2);
0214 isRHBehindECAL &= dRho2 <= dRho2Max_HCALClus_SamePhi_;
0215 if (isRHBehindECAL) {
0216 returnValue.samedPhi_.x_ += rhX * rhE;
0217 returnValue.samedPhi_.y_ += rhY * rhE;
0218 returnValue.samedPhi_.z_ += rhZ * rhE;
0219 returnValue.samedPhi_.e_ += rhE;
0220 returnValue.samedPhi_.nHits_++;
0221 }
0222 }
0223 }
0224
0225
0226 if (!isRHBehindECAL) {
0227 double dR2 = reco::deltaR2(phoSCEta, phoSCPhi, rhEta, rhPhi);
0228 if (dR2 < dR2Max_HCALClus_SamePhi_) {
0229 returnValue.samedR_.x_ += rhX * rhE;
0230 returnValue.samedR_.y_ += rhY * rhE;
0231 returnValue.samedR_.z_ += rhZ * rhE;
0232 returnValue.samedR_.e_ += rhE;
0233 returnValue.samedR_.nHits_++;
0234 }
0235 }
0236 }
0237
0238 if (returnValue.samedPhi_.nHits_ > 0) {
0239 returnValue.samedPhi_.x_ = returnValue.samedPhi_.x_ / returnValue.samedPhi_.e_;
0240 returnValue.samedPhi_.y_ = returnValue.samedPhi_.y_ / returnValue.samedPhi_.e_;
0241 returnValue.samedPhi_.z_ = returnValue.samedPhi_.z_ / returnValue.samedPhi_.e_;
0242 }
0243
0244 if (returnValue.samedR_.nHits_ > 0) {
0245 returnValue.samedR_.x_ = returnValue.samedR_.x_ / returnValue.samedR_.e_;
0246 returnValue.samedR_.y_ = returnValue.samedR_.y_ / returnValue.samedR_.e_;
0247 returnValue.samedR_.z_ = returnValue.samedR_.z_ / returnValue.samedR_.e_;
0248 }
0249 return returnValue;
0250 }
0251
0252 PhotonMVABasedHaloTagger::PreshowerHyp PhotonMVABasedHaloTagger::calmatchedESCoordForBothHypothesis(
0253 const CaloGeometry& geo,
0254 const reco::Photon* pho,
0255 const EcalRecHitCollection& ESRecHits,
0256 const EcalClus& ecalClus) const {
0257 PreshowerHyp returnValue;
0258
0259 double phoSCEta = pho->superCluster()->eta();
0260 double phoSCPhi = pho->superCluster()->phi();
0261
0262 double tmpDiffdRho = 999;
0263 double matchX_samephi = -999;
0264 double matchY_samephi = -999;
0265 bool foundESRH_samephi = false;
0266
0267 double tmpDiffdRho_samedR = 999;
0268 double matchX_samedR = -999;
0269 double matchY_samedR = -999;
0270 bool foundESRH_samedR = false;
0271
0272
0273 double tan_theta = 1. / sinh(phoSCEta);
0274
0275 double cos_phi = cos(phoSCPhi);
0276 double sin_phi = sin(phoSCPhi);
0277
0278 for (const auto& esrechit : ESRecHits) {
0279 const GlobalPoint& rechitPoint = geo.getPosition(esrechit.id());
0280
0281 double rhEta = rechitPoint.eta();
0282 double rhX = rechitPoint.x();
0283 double rhY = rechitPoint.y();
0284 double rhZ = rechitPoint.z();
0285 double rhE = esrechit.energy();
0286
0287 if (phoSCEta * rhEta < 0)
0288 continue;
0289
0290 if (rhE < noiseThrES_)
0291 continue;
0292
0293
0294
0295
0296
0297
0298 double dRho2 = pow(rhX - ecalClus.x_, 2) + pow(rhY - ecalClus.y_, 2);
0299
0300 if (dRho2 < tmpDiffdRho && dRho2 < dRho2Max_ESClus_) {
0301 tmpDiffdRho = dRho2;
0302 matchX_samephi = rhX;
0303 matchY_samephi = rhY;
0304 foundESRH_samephi = true;
0305 }
0306
0307
0308
0309 double exp_ESRho = rhZ * tan_theta;
0310 double exp_ESX = cos_phi * exp_ESRho;
0311 double exp_ESY = sin_phi * exp_ESRho;
0312
0313 double dRho_samedR2 = pow(rhX - exp_ESX, 2) + pow(rhY - exp_ESY, 2);
0314
0315 if (dRho_samedR2 < tmpDiffdRho_samedR) {
0316 tmpDiffdRho_samedR = dRho_samedR2;
0317 matchX_samedR = rhX;
0318 matchY_samedR = rhY;
0319 foundESRH_samedR = true;
0320 }
0321
0322 }
0323
0324
0325
0326
0327 for (const auto& esrechit : ESRecHits) {
0328 const GlobalPoint& rechitPoint = geo.getPosition(esrechit.id());
0329
0330 double rhEta = rechitPoint.eta();
0331 double rhX = rechitPoint.x();
0332 double rhY = rechitPoint.y();
0333 double rhZ = rechitPoint.z();
0334 double rhE = esrechit.energy();
0335
0336 if (phoSCEta * rhEta < 0)
0337 continue;
0338 if (rhE < noiseThrES_)
0339 continue;
0340
0341
0342 bool isRHBehindECAL = foundESRH_samephi;
0343 if (isRHBehindECAL) {
0344 double dX_samephi = std::abs(matchX_samephi - rhX);
0345 double dY_samephi = std::abs(matchY_samephi - rhY);
0346 isRHBehindECAL &= (dX_samephi < dXY_ESClus_SamePhi_ && dY_samephi < dXY_ESClus_SamePhi_);
0347 if (isRHBehindECAL) {
0348 returnValue.samedPhi_.x_ += rhX * rhE;
0349 returnValue.samedPhi_.y_ += rhY * rhE;
0350 returnValue.samedPhi_.z_ += rhZ * rhE;
0351 returnValue.samedPhi_.e_ += rhE;
0352 returnValue.samedPhi_.nHits_++;
0353 }
0354 }
0355
0356
0357 if (!isRHBehindECAL && foundESRH_samedR) {
0358 double dX_samedR = std::abs(matchX_samedR - rhX);
0359 double dY_samedR = std::abs(matchY_samedR - rhY);
0360
0361 if (dX_samedR < dXY_ESClus_SamedR_ && dY_samedR < dXY_ESClus_SamedR_) {
0362 returnValue.samedR_.x_ += rhX * rhE;
0363 returnValue.samedR_.y_ += rhY * rhE;
0364 returnValue.samedR_.z_ += rhZ * rhE;
0365 returnValue.samedR_.e_ += rhE;
0366 returnValue.samedR_.nHits_++;
0367 }
0368 }
0369 }
0370
0371 if (returnValue.samedPhi_.nHits_ > 0) {
0372 returnValue.samedPhi_.x_ = returnValue.samedPhi_.x_ / returnValue.samedPhi_.e_;
0373 returnValue.samedPhi_.y_ = returnValue.samedPhi_.y_ / returnValue.samedPhi_.e_;
0374 returnValue.samedPhi_.z_ = returnValue.samedPhi_.z_ / returnValue.samedPhi_.e_;
0375 }
0376
0377 if (returnValue.samedR_.nHits_ > 0) {
0378 returnValue.samedR_.x_ = returnValue.samedR_.x_ / returnValue.samedR_.e_;
0379 returnValue.samedR_.y_ = returnValue.samedR_.y_ / returnValue.samedR_.e_;
0380 returnValue.samedR_.z_ = returnValue.samedR_.z_ / returnValue.samedR_.e_;
0381 }
0382 return returnValue;
0383 }
0384
0385 double PhotonMVABasedHaloTagger::calAngleBetweenEEAndSubDet(CalClus const& subdetClus, EcalClus const& ecalClus) const {
0386
0387
0388 double angle = -999;
0389
0390 if (ecalClus.nHits_ > 0 && subdetClus.nHits_ > 0) {
0391 double dR = sqrt(pow(subdetClus.x_ - ecalClus.x_, 2) + pow(subdetClus.y_ - ecalClus.y_, 2) +
0392 pow(subdetClus.z_ - ecalClus.z_, 2));
0393
0394 double cosTheta = std::abs(subdetClus.z_ - ecalClus.z_) / dR;
0395
0396 angle = acos(cosTheta);
0397 }
0398
0399 return angle;
0400 }