PhotonIdentification/src/PhotonMVABasedHaloTagger.cc

0001 /** \class PhotonMVABasedHaloTagger
0002  *   \author Shilpi Jain (University of Minnesota)
0003  *  * Links to the presentation:
0004  1. ECAL DPG: https://indico.cern.ch/event/991261/contributions/4283096/attachments/2219229/3757719/beamHalo_31march_v1.pdf
0005  2. JetMET POG: https://indico.cern.ch/event/1027614/contributions/4314949/attachments/2224472/3767396/beamHalo_12April.pdf
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()),
0020       ecalPFRechitThresholdsToken_(iC.esConsumes()),
0021       ecalClusterToolsESGetTokens_(std::move(iC)) {
0022   rhoLabel_ = iC.consumes<double>(conf.getParameter<edm::InputTag>("rhoLabel"));
0023
0024   EBecalCollection_ = iC.consumes<EcalRecHitCollection>(conf.getParameter<edm::InputTag>("barrelEcalRecHitCollection"));
0025   EEecalCollection_ = iC.consumes<EcalRecHitCollection>(conf.getParameter<edm::InputTag>("endcapEcalRecHitCollection"));
0026   ESCollection_ = iC.consumes<EcalRecHitCollection>(conf.getParameter<edm::InputTag>("esRecHitCollection"));
0027   HBHERecHitsCollection_ = iC.consumes<HBHERecHitCollection>(conf.getParameter<edm::InputTag>("HBHERecHitsCollection"));
0028   recHitEThresholdHB_ = conf.getParameter<EgammaHcalIsolation::arrayHB>("recHitEThresholdHB");
0029   recHitEThresholdHE_ = conf.getParameter<EgammaHcalIsolation::arrayHE>("recHitEThresholdHE");
0030
0031   noiseThrES_ = conf.getParameter<double>("noiseThrES");
0032 }
0033
0034 double PhotonMVABasedHaloTagger::calculateMVA(const reco::Photon* pho,
0035                                               const GBRForest* gbr_,
0036                                               const edm::Event& iEvent,
0037                                               const edm::EventSetup& es) {
0038   bool isEB = pho->isEB();
0039
0040   if (isEB)
0041     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
0042
0043   //rho handle
0044   double rho_ = iEvent.get(rhoLabel_);
0045
0046   // Get all the RecHits
0047   const auto& ecalRecHitsBarrel = iEvent.get(EBecalCollection_);
0048   const auto& ecalRecHitsEndcap = iEvent.get(EEecalCollection_);
0049   const auto& esRecHits = iEvent.get(ESCollection_);
0050   const auto& hbheRecHits = iEvent.get(HBHERecHitsCollection_);
0051
0052   //gets geometry
0053   pG_ = es.getHandle(geometryToken_);
0054   const CaloGeometry* geo = pG_.product();
0055
0056   ///ECAL PF rechit thresholds
0057   auto const& thresholds = es.getData(ecalPFRechitThresholdsToken_);
0058
0059   noZS::EcalClusterLazyTools lazyToolnoZS(
0060       iEvent, ecalClusterToolsESGetTokens_.get(es), EBecalCollection_, EEecalCollection_);
0061
0062   ///calculate the energy weighted X, Y and Z position of the photon cluster
0063   if (isEB)
0064     calphoClusCoordinECAL(geo, pho, &thresholds, ecalRecHitsBarrel);
0065   else
0066     calphoClusCoordinECAL(geo, pho, &thresholds, ecalRecHitsEndcap);
0067
0068   ///calculate the HBHE cluster position hypothesis
0069   calmatchedHBHECoordForBothHypothesis(geo, pho, hbheRecHits);
0070   calmatchedESCoordForBothHypothesis(geo, pho, esRecHits);
0071
0072   ///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
0073
0074   double angle_EE_HE_samedPhi = calAngleBetweenEEAndSubDet(
0075       hcalClusNhits_samedPhi_,
0076       hcalClusX_samedPhi_,
0077       hcalClusY_samedPhi_,
0078       hcalClusZ_samedPhi_);  //essentially caculates the angle and energy variables in the two hypothesis between EE and HE
0079
0080   double angle_EE_HE_samedR =
0081       calAngleBetweenEEAndSubDet(hcalClusNhits_samedR_, hcalClusX_samedR_, hcalClusY_samedR_, hcalClusZ_samedR_);
0082
0083   double angle_EE_ES_samedPhi = calAngleBetweenEEAndSubDet(
0084       preshowerNhits_samedPhi_, preshowerX_samedPhi_, preshowerY_samedPhi_, preshowerZ_samedPhi_);
0085
0086   double angle_EE_ES_samedR =
0087       calAngleBetweenEEAndSubDet(preshowerNhits_samedR_, preshowerX_samedR_, preshowerY_samedR_, preshowerZ_samedR_);
0088
0089   ////set all the above calculated variables as input to the MVA
0090
0091   const auto& vCov = lazyToolnoZS.localCovariances(*(pho->superCluster()->seed()));
0092   double spp = (std::isnan(vCov[2]) ? 0. : sqrt(vCov[2]));
0093
0094   ///https://cmssdt.cern.ch/lxr/source/RecoEgamma/ElectronIdentification/src/ElectronMVAEstimator.cc
0095
0096   float vars[15];
0097
0098   vars[0] = preshowerE_samedPhi_;
0099   vars[1] = hcalClusE_samedPhi_;
0100   vars[2] = preshowerE_samedR_;
0101   vars[3] = hcalClusE_samedR_;
0102   vars[4] = pho->full5x5_r9();
0103   vars[5] = pho->superCluster()->etaWidth();
0104   vars[6] = pho->superCluster()->phiWidth();
0105   vars[7] = pho->full5x5_sigmaIetaIeta();
0106   vars[8] = spp;
0107   vars[9] = angle_EE_ES_samedR;
0108   vars[10] = angle_EE_HE_samedR;
0109   vars[11] = angle_EE_ES_samedPhi;
0110   vars[12] = angle_EE_HE_samedPhi;
0111   vars[13] = (pho->superCluster()->preshowerEnergyPlane1() + pho->superCluster()->preshowerEnergyPlane2()) /
0112              pho->superCluster()->rawEnergy();
0113   vars[14] = rho_;
0114
0115   double BHmva = gbr_->GetGradBoostClassifier(vars);
0116   return BHmva;
0117 }
0118
0119 void PhotonMVABasedHaloTagger::calphoClusCoordinECAL(const CaloGeometry* geo,
0120                                                      const reco::Photon* pho,
0121                                                      const EcalPFRecHitThresholds* thresholds,
0122                                                      const EcalRecHitCollection& ecalRecHits) {
0123   ecalClusX_ = 0;
0124   ecalClusY_ = 0;
0125   ecalClusZ_ = 0;
0126   ecalClusNhits_ = 0;
0127   ecalClusE_ = 0;
0128
0129   double phoSCEta = pho->superCluster()->eta();
0130   double phoSCPhi = pho->superCluster()->phi();
0131
0132   for (const auto& ecalrechit : ecalRecHits) {
0133     auto const det = ecalrechit.id();
0134     double rhE = ecalrechit.energy();
0135     const GlobalPoint& rechitPoint = geo->getPosition(det);
0136
0137     double rhEta = rechitPoint.eta();
0138     double rhPhi = rechitPoint.phi();
0139     double rhX = rechitPoint.x();
0140     double rhY = rechitPoint.y();
0141     double rhZ = rechitPoint.z();
0142
0143     if (thresholds == nullptr) {
0144       throw cms::Exception("EmptyPFRechHitThresCollection")
0145           << "In PhotonMVABasedHaloTagger::calphoClusCoordinECAL, EcalPFRecHitThresholds cannot be = nulptr";
0146     }
0147
0148     float rhThres = (*thresholds)[det];
0149
0150     if (rhE <= rhThres)
0151       continue;
0152
0153     if (phoSCEta * rhEta < 0)
0154       continue;
0155
0156     double dR2 = reco::deltaR2(rhEta, rhPhi, phoSCEta, phoSCPhi);
0157
0158     if (dR2 < dr2Max_ECALClus_) {
0159       ecalClusX_ += rhX * rhE;
0160       ecalClusY_ += rhY * rhE;
0161       ecalClusZ_ += rhZ * rhE;
0162       ecalClusE_ += rhE;
0163       ecalClusNhits_++;
0164     }
0165   }  //for(int ih=0; ih<nMatchedEErh[ipho]; ih++)
0166
0167   if (ecalClusNhits_ > 0) {  //should always be > 0 for an EM cluster
0168     ecalClusX_ = ecalClusX_ / ecalClusE_;
0169     ecalClusY_ = ecalClusY_ / ecalClusE_;
0170     ecalClusZ_ = ecalClusZ_ / ecalClusE_;
0171   }  //if(ecalClusNhits_>0)
0172 }
0173
0174 void PhotonMVABasedHaloTagger::calmatchedHBHECoordForBothHypothesis(const CaloGeometry* geo,
0175                                                                     const reco::Photon* pho,
0176                                                                     const HBHERecHitCollection& HBHERecHits) {
0177   hcalClusX_samedPhi_ = 0;
0178   hcalClusY_samedPhi_ = 0;
0179   hcalClusZ_samedPhi_ = 0;
0180
0181   hcalClusNhits_samedPhi_ = 0;
0182   hcalClusE_samedPhi_ = 0;
0183
0184   hcalClusX_samedR_ = 0;
0185   hcalClusY_samedR_ = 0;
0186   hcalClusZ_samedR_ = 0;
0187   hcalClusNhits_samedR_ = 0;
0188   hcalClusE_samedR_ = 0;
0189
0190   double phoSCEta = pho->superCluster()->eta();
0191   double phoSCPhi = pho->superCluster()->phi();
0192
0193   // Loop over HBHERecHit's
0194   for (const auto& hbherechit : HBHERecHits) {
0195     HcalDetId det = hbherechit.id();
0196     const GlobalPoint& rechitPoint = geo->getPosition(det);
0197
0198     double rhEta = rechitPoint.eta();
0199     double rhPhi = rechitPoint.phi();
0200     double rhX = rechitPoint.x();
0201     double rhY = rechitPoint.y();
0202     double rhZ = rechitPoint.z();
0203     double rhE = hbherechit.energy();
0204
0205     int depth = det.depth();
0206
0207     if ((det.subdet() == HcalBarrel and (depth < 1 or depth > int(recHitEThresholdHB_.size()))) or
0208         (det.subdet() == HcalEndcap and (depth < 1 or depth > int(recHitEThresholdHE_.size())))) {
0209       edm::LogWarning("PhotonMVABasedHaloTagger")
0210           << " hit in subdet " << det.subdet() << " has an unaccounted for depth of " << depth
0211           << "!! Leaving this hit!!";
0212       continue;
0213     }
0214
0215     const bool goodHBe = det.subdet() == HcalBarrel and rhE > recHitEThresholdHB_[depth - 1];
0216     const bool goodHEe = det.subdet() == HcalEndcap and rhE > recHitEThresholdHE_[depth - 1];
0217     if (!(goodHBe or goodHEe))
0218       continue;
0219
0220     if (phoSCEta * rhEta < 0)
0221       continue;  ///Should be on the same side of Z
0222
0223     double dPhi = deltaPhi(phoSCPhi, rhPhi);
0224
0225     ///only valid for the EE; this is 26 cm; hit within 3x3 of HCAL centered at the EECAL xtal
0226     bool isRHBehindECAL = std::abs(dPhi) < dPhiMax_HCALClus_SamePhi_;
0227     if (isRHBehindECAL) {
0228       double rho2 = pow(rhX, 2) + pow(rhY, 2);
0229       isRHBehindECAL &= (rho2 >= rho2Min_ECALpos_ && rho2 <= rho2Max_ECALpos_);
0230       if (isRHBehindECAL) {
0231         double dRho2 = pow(rhX - ecalClusX_, 2) + pow(rhY - ecalClusY_, 2);
0232         isRHBehindECAL &= dRho2 <= dRho2Max_HCALClus_SamePhi_;
0233         if (isRHBehindECAL) {
0234           hcalClusX_samedPhi_ += rhX * rhE;
0235           hcalClusY_samedPhi_ += rhY * rhE;
0236           hcalClusZ_samedPhi_ += rhZ * rhE;
0237           hcalClusE_samedPhi_ += rhE;
0238           hcalClusNhits_samedPhi_++;
0239         }
0240       }
0241     }  //if(rho>=31 && rho<=172)
0242
0243     ///dont use hits which are just behind the ECAL in the same phi region
0244     if (!isRHBehindECAL) {
0245       double dR2 = reco::deltaR2(phoSCEta, phoSCPhi, rhEta, rhPhi);
0246       if (dR2 < dR2Max_HCALClus_SamePhi_) {
0247         hcalClusX_samedR_ += rhX * rhE;
0248         hcalClusY_samedR_ += rhY * rhE;
0249         hcalClusZ_samedR_ += rhZ * rhE;
0250         hcalClusE_samedR_ += rhE;
0251         hcalClusNhits_samedR_++;
0252       }
0253     }
0254   }  //for(int ih=0; ih<nMatchedHBHErh[ipho]; ih++)
0255
0256   if (hcalClusNhits_samedPhi_ > 0) {
0257     hcalClusX_samedPhi_ = hcalClusX_samedPhi_ / hcalClusE_samedPhi_;
0258     hcalClusY_samedPhi_ = hcalClusY_samedPhi_ / hcalClusE_samedPhi_;
0259     hcalClusZ_samedPhi_ = hcalClusZ_samedPhi_ / hcalClusE_samedPhi_;
0260   }  //if(hcalClusNhits_samedPhi_>0)
0261
0262   if (hcalClusNhits_samedR_ > 0) {
0263     hcalClusX_samedR_ = hcalClusX_samedR_ / hcalClusE_samedR_;
0264     hcalClusY_samedR_ = hcalClusY_samedR_ / hcalClusE_samedR_;
0265     hcalClusZ_samedR_ = hcalClusZ_samedR_ / hcalClusE_samedR_;
0266   }  //if(hcalClusNhits_samedR_>0)
0267 }
0268
0269 void PhotonMVABasedHaloTagger::calmatchedESCoordForBothHypothesis(const CaloGeometry* geo,
0270                                                                   const reco::Photon* pho,
0271                                                                   const EcalRecHitCollection& ESRecHits) {
0272   preshowerX_samedPhi_ = 0;
0273   preshowerY_samedPhi_ = 0;
0274   preshowerZ_samedPhi_ = 0;
0275   preshowerNhits_samedPhi_ = 0;
0276   preshowerE_samedPhi_ = 0;
0277
0278   preshowerX_samedR_ = 0;
0279   preshowerY_samedR_ = 0;
0280   preshowerZ_samedR_ = 0;
0281   preshowerNhits_samedR_ = 0;
0282   preshowerE_samedR_ = 0;
0283
0284   double phoSCEta = pho->superCluster()->eta();
0285   double phoSCPhi = pho->superCluster()->phi();
0286
0287   double tmpDiffdRho = 999;
0288   double matchX_samephi = -999;
0289   double matchY_samephi = -999;
0290   bool foundESRH_samephi = false;
0291
0292   double tmpDiffdRho_samedR = 999;
0293   double matchX_samedR = -999;
0294   double matchY_samedR = -999;
0295   bool foundESRH_samedR = false;
0296
0297   ///get theta and phi of the coordinates of photon
0298   double tan_theta = 1. / sinh(phoSCEta);
0299
0300   double cos_phi = cos(phoSCPhi);
0301   double sin_phi = sin(phoSCPhi);
0302
0303   for (const auto& esrechit : ESRecHits) {
0304     const GlobalPoint& rechitPoint = geo->getPosition(esrechit.id());
0305
0306     double rhEta = rechitPoint.eta();
0307     double rhX = rechitPoint.x();
0308     double rhY = rechitPoint.y();
0309     double rhZ = rechitPoint.z();
0310     double rhE = esrechit.energy();
0311
0312     if (phoSCEta * rhEta < 0)
0313       continue;
0314
0315     if (rhE < noiseThrES_)
0316       continue;
0317
0318     ////try to include RH according to the strips, 11 in X and 11 in Y
0319     /////////First calculate RH nearest in phi and eta to that of the photon SC
0320
0321     //////same phi ----> the X and Y should be similar
0322     ////i.e. hit is required to be within the ----> seems better match with the data compared to 2.47
0323     double dRho2 = pow(rhX - ecalClusX_, 2) + pow(rhY - ecalClusY_, 2);
0324
0325     if (dRho2 < tmpDiffdRho && dRho2 < dRho2Max_ESClus_) {
0326       tmpDiffdRho = dRho2;
0327       matchX_samephi = rhX;
0328       matchY_samephi = rhY;
0329       foundESRH_samephi = true;
0330     }
0331
0332     ////////same eta
0333     ///calculate the expected x and y at the position of hte rechit
0334     double exp_ESRho = rhZ * tan_theta;
0335     double exp_ESX = cos_phi * exp_ESRho;
0336     double exp_ESY = sin_phi * exp_ESRho;
0337
0338     double dRho_samedR2 = pow(rhX - exp_ESX, 2) + pow(rhY - exp_ESY, 2);
0339
0340     if (dRho_samedR2 < tmpDiffdRho_samedR) {
0341       tmpDiffdRho_samedR = dRho_samedR2;
0342       matchX_samedR = rhX;
0343       matchY_samedR = rhY;
0344       foundESRH_samedR = true;
0345     }
0346
0347   }  ///  for (const auto& esrechit : ESRecHits)
0348
0349   ////Now calculate the sum in +/- 5 strips in X and y around the matched RH
0350   //+/5  strips mean = 5*~2mm = +/-10 mm = 1 cm
0351
0352   for (const auto& esrechit : ESRecHits) {
0353     const GlobalPoint& rechitPoint = geo->getPosition(esrechit.id());
0354
0355     double rhEta = rechitPoint.eta();
0356     double rhX = rechitPoint.x();
0357     double rhY = rechitPoint.y();
0358     double rhZ = rechitPoint.z();
0359     double rhE = esrechit.energy();
0360
0361     if (phoSCEta * rhEta < 0)
0362       continue;
0363     if (rhE < noiseThrES_)
0364       continue;
0365
0366     ///same phi
0367     bool isRHBehindECAL = foundESRH_samephi;
0368     if (isRHBehindECAL) {
0369       double dX_samephi = std::abs(matchX_samephi - rhX);
0370       double dY_samephi = std::abs(matchY_samephi - rhY);
0371       isRHBehindECAL &= (dX_samephi < dXY_ESClus_SamePhi_ && dY_samephi < dXY_ESClus_SamePhi_);
0372       if (isRHBehindECAL) {
0373         preshowerX_samedPhi_ += rhX * rhE;
0374         preshowerY_samedPhi_ += rhY * rhE;
0375         preshowerZ_samedPhi_ += rhZ * rhE;
0376         preshowerE_samedPhi_ += rhE;
0377         preshowerNhits_samedPhi_++;
0378       }
0379     }
0380
0381     ///same dR
0382     if (!isRHBehindECAL && foundESRH_samedR) {
0383       double dX_samedR = std::abs(matchX_samedR - rhX);
0384       double dY_samedR = std::abs(matchY_samedR - rhY);
0385
0386       if (dX_samedR < dXY_ESClus_SamedR_ && dY_samedR < dXY_ESClus_SamedR_) {
0387         preshowerX_samedR_ += rhX * rhE;
0388         preshowerY_samedR_ += rhY * rhE;
0389         preshowerZ_samedR_ += rhZ * rhE;
0390         preshowerE_samedR_ += rhE;
0391         preshowerNhits_samedR_++;
0392       }
0393     }
0394   }  ///for(int ih=0; ih<nMatchedESrh[ipho]; ih++)
0395
0396   if (preshowerNhits_samedPhi_ > 0) {
0397     preshowerX_samedPhi_ = preshowerX_samedPhi_ / preshowerE_samedPhi_;
0398     preshowerY_samedPhi_ = preshowerY_samedPhi_ / preshowerE_samedPhi_;
0399     preshowerZ_samedPhi_ = preshowerZ_samedPhi_ / preshowerE_samedPhi_;
0400   }  //if(preshowerNhits_samedPhi_>0)
0401
0402   if (preshowerNhits_samedR_ > 0) {
0403     preshowerX_samedR_ = preshowerX_samedR_ / preshowerE_samedR_;
0404     preshowerY_samedR_ = preshowerY_samedR_ / preshowerE_samedR_;
0405     preshowerZ_samedR_ = preshowerZ_samedR_ / preshowerE_samedR_;
0406   }  //if(preshowerNhits_samedR_>0)
0407 }
0408
0409 double PhotonMVABasedHaloTagger::calAngleBetweenEEAndSubDet(int nhits,
0410                                                             double subdetClusX,
0411                                                             double subdetClusY,
0412                                                             double subdetClusZ) {
0413   ////get the angle of the line joining the ECAL cluster and the subdetector wrt Z axis for any hypothesis
0414
0415   double angle = -999;
0416
0417   if (ecalClusNhits_ > 0 && nhits > 0) {
0418     double dR =
0419         sqrt(pow(subdetClusX - ecalClusX_, 2) + pow(subdetClusY - ecalClusY_, 2) + pow(subdetClusZ - ecalClusZ_, 2));
0420
0421     double cosTheta = std::abs(subdetClusZ - ecalClusZ_) / dR;
0422
0423     angle = acos(cosTheta);
0424   }
0425
0426   return angle;
0427 }