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()), 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;  /// 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
0040
0041   //rho handle
0042   double rho_ = iEvent.get(rhoLabel_);
0043
0044   // Get all the RecHits
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   //gets geometry
0051   const CaloGeometry& geo = es.getData(geometryToken_);
0052
0053   ///ECAL PF rechit thresholds
0054   auto const& thresholds = es.getData(ecalPFRechitThresholdsToken_);
0055
0056   noZS::EcalClusterLazyTools lazyToolnoZS(
0057       iEvent, ecalClusterToolsESGetTokens_.get(es), EBecalCollection_, EEecalCollection_);
0058
0059   ///calculate the energy weighted X, Y and Z position of the photon cluster
0060   EcalClus ecalClus;
0061   if (isEB)
0062     ecalClus = calphoClusCoordinECAL(geo, pho, &thresholds, ecalRecHitsBarrel);
0063   else
0064     ecalClus = calphoClusCoordinECAL(geo, pho, &thresholds, ecalRecHitsEndcap);
0065
0066   ///calculate the HBHE cluster position hypothesis
0067   auto hcalClus = calmatchedHBHECoordForBothHypothesis(geo, pho, hbheRecHits, ecalClus);
0068   auto preshower = calmatchedESCoordForBothHypothesis(geo, pho, esRecHits, ecalClus);
0069
0070   ///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
0071
0072   double angle_EE_HE_samedPhi = calAngleBetweenEEAndSubDet(
0073       hcalClus.samedPhi_,
0074       ecalClus);  //essentially caculates the angle and energy variables in the two hypothesis between EE and HE
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   ////set all the above calculated variables as input to the MVA
0083
0084   const auto& vCov = lazyToolnoZS.localCovariances(*(pho->superCluster()->seed()));
0085   double spp = (std::isnan(vCov[2]) ? 0. : sqrt(vCov[2]));
0086
0087   ///https://cmssdt.cern.ch/lxr/source/RecoEgamma/ElectronIdentification/src/ElectronMVAEstimator.cc
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   }  //for(int ih=0; ih<nMatchedEErh[ipho]; ih++)
0155
0156   if (ecalClus.nHits_ > 0) {  //should always be > 0 for an EM cluster
0157     ecalClus.x_ = ecalClus.x_ / ecalClus.e_;
0158     ecalClus.y_ = ecalClus.y_ / ecalClus.e_;
0159     ecalClus.z_ = ecalClus.z_ / ecalClus.e_;
0160   }  //if(ecalClus.nHits_>0)
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   // Loop over HBHERecHit's
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;  ///Should be on the same side of Z
0204
0205     double dPhi = deltaPhi(phoSCPhi, rhPhi);
0206
0207     ///only valid for the EE; this is 26 cm; hit within 3x3 of HCAL centered at the EECAL xtal
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     }  //if(rho>=31 && rho<=172)
0224
0225     ///dont use hits which are just behind the ECAL in the same phi region
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   }  //for(int ih=0; ih<nMatchedHBHErh[ipho]; ih++)
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   }  //if(returnValue.samedPhi_.Nhits_>0)
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   }  //if(returnValue.samedR_.nHits_>0)
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   ///get theta and phi of the coordinates of photon
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     ////try to include RH according to the strips, 11 in X and 11 in Y
0294     /////////First calculate RH nearest in phi and eta to that of the photon SC
0295
0296     //////same phi ----> the X and Y should be similar
0297     ////i.e. hit is required to be within the ----> seems better match with the data compared to 2.47
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     ////////same eta
0308     ///calculate the expected x and y at the position of hte rechit
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   }  ///  for (const auto& esrechit : ESRecHits)
0323
0324   ////Now calculate the sum in +/- 5 strips in X and y around the matched RH
0325   //+/5  strips mean = 5*~2mm = +/-10 mm = 1 cm
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     ///same phi
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     ///same dR
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   }  ///for(int ih=0; ih<nMatchedESrh[ipho]; ih++)
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   }  //if(preshowerSamedPhi_.nHits_>0)
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   }  //if(preshowerSamedR_.nHits_>0)
0382   return returnValue;
0383 }
0384
0385 double PhotonMVABasedHaloTagger::calAngleBetweenEEAndSubDet(CalClus const& subdetClus, EcalClus const& ecalClus) const {
0386   ////get the angle of the line joining the ECAL cluster and the subdetector wrt Z axis for any hypothesis
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 }