EgammaIsolationAlgos/src/EgammaRecHitIsolation.cc

0001 //*****************************************************************************
0002 // File:      EgammaRecHitIsolation.cc
0003 // ----------------------------------------------------------------------------
0004 // OrigAuth:  Matthias Mozer, hacked by Sam Harper (ie the ugly stuff is mine)
0005 // Institute: IIHE-VUB, RAL
0006 //=============================================================================
0007 //*****************************************************************************
0008 //C++ includes
0009 #include <vector>
0010 #include <functional>
0011
0012 //ROOT includes
0013 #include <Math/VectorUtil.h>
0014
0015 //CMSSW includes
0016 #include "RecoEgamma/EgammaIsolationAlgos/interface/EgammaRecHitIsolation.h"
0017 #include "DataFormats/DetId/interface/DetId.h"
0018 #include "DataFormats/EcalDetId/interface/EcalSubdetector.h"
0019 #include "DataFormats/GeometryVector/interface/GlobalPoint.h"
0020 #include "DataFormats/GeometryVector/interface/GlobalVector.h"
0021 #include "DataFormats/EgammaReco/interface/BasicCluster.h"
0022 #include "DataFormats/EgammaReco/interface/BasicClusterFwd.h"
0023 #include "DataFormats/CaloRecHit/interface/CaloCluster.h"
0024 #include "DataFormats/CaloRecHit/interface/CaloClusterFwd.h"
0025 #include "DataFormats/EgammaReco/interface/SuperCluster.h"
0026 #include "DataFormats/EgammaReco/interface/SuperClusterFwd.h"
0027 #include "DataFormats/RecoCandidate/interface/RecoCandidate.h"
0028 #include "DataFormats/Math/interface/deltaPhi.h"
0029
0030 using namespace std;
0031
0032 EgammaRecHitIsolation::EgammaRecHitIsolation(double extRadius,
0033                                              double intRadius,
0034                                              double etaSlice,
0035                                              double etLow,
0036                                              double eLow,
0037                                              edm::ESHandle<CaloGeometry> theCaloGeom,
0038                                              const EcalRecHitCollection& caloHits,
0039                                              const EcalSeverityLevelAlgo* sl,
0040                                              DetId::Detector detector)
0041     :  // not used anymore, kept for compatibility
0042       extRadius_(extRadius),
0043       intRadius_(intRadius),
0044       etaSlice_(etaSlice),
0045       etLow_(etLow),
0046       eLow_(eLow),
0047       theCaloGeom_(theCaloGeom),
0048       caloHits_(caloHits),
0049       sevLevel_(sl),
0050       useNumCrystals_(false),
0051       vetoClustered_(false),
0052       ecalBarHits_(nullptr),
0053       //chStatus_(0),
0054       severitiesexcl_(0),
0055       //severityRecHitThreshold_(0),
0056       //spId_(EcalSeverityLevelAlgo::kSwissCross),
0057       //spIdThreshold_(0),
0058       flags_(0) {
0059   //set up the geometry and selector
0060   const CaloGeometry* caloGeom = theCaloGeom_.product();
0061   subdet_[0] = caloGeom->getSubdetectorGeometry(DetId::Ecal, EcalBarrel);
0062   subdet_[1] = caloGeom->getSubdetectorGeometry(DetId::Ecal, EcalEndcap);
0063 }
0064
0065 EgammaRecHitIsolation::~EgammaRecHitIsolation() {}
0066
0067 double EgammaRecHitIsolation::getSum_(const reco::Candidate* emObject,
0068                                       bool returnEt,
0069                                       const EcalPFRecHitThresholds* thresholds) const {
0070   double energySum = 0.;
0071   if (!caloHits_.empty()) {
0072     //Take the SC position
0073     reco::SuperClusterRef sc = emObject->get<reco::SuperClusterRef>();
0074     math::XYZPoint const& theCaloPosition = sc.get()->position();
0075     GlobalPoint pclu(theCaloPosition.x(), theCaloPosition.y(), theCaloPosition.z());
0076     float etaclus = pclu.eta();
0077     float phiclus = pclu.phi();
0078     float r2 = intRadius_ * intRadius_;
0079
0080     std::vector<std::pair<DetId, float> >::const_iterator rhIt;
0081
0082     for (int subdetnr = 0; subdetnr <= 1; subdetnr++) {  // look in barrel and endcap
0083       if (nullptr == subdet_[subdetnr])
0084         continue;
0085
0086       CaloSubdetectorGeometry::DetIdSet chosen =
0087           subdet_[subdetnr]->getCells(pclu, extRadius_);  // select cells around cluster
0088       EcalRecHitCollection::const_iterator j = caloHits_.end();
0089
0090       for (CaloSubdetectorGeometry::DetIdSet::const_iterator i = chosen.begin(); i != chosen.end();
0091            ++i) {                    //loop selected cells
0092         j = caloHits_.find(*i);      // find selected cell among rechits
0093         if (j != caloHits_.end()) {  // add rechit only if available
0094           auto cell = theCaloGeom_->getGeometry(*i);
0095           float eta = cell->etaPos();
0096           float phi = cell->phiPos();
0097           float etaDiff = eta - etaclus;
0098           float phiDiff = reco::deltaPhi(phi, phiclus);
0099           float energy = j->energy();
0100
0101           float rhThres = (thresholds != nullptr) ? (*thresholds)[j->detid()] : 0.f;
0102           if (energy <= rhThres)
0103             continue;
0104
0105           if (useNumCrystals_) {
0106             if (fabs(etaclus) < 1.479) {  // Barrel num crystals, crystal width = 0.0174
0107               if (fabs(etaDiff) < 0.0174 * etaSlice_)
0108                 continue;
0109               //if (sqrt(etaDiff*etaDiff + phiDiff*phiDiff) < 0.0174*intRadius_)
0110               //continue;
0111               if ((etaDiff * etaDiff + phiDiff * phiDiff) < 0.00030276 * r2)
0112                 continue;
0113             } else {  // Endcap num crystals, crystal width = 0.00864*fabs(sinh(eta))
0114               if (fabs(etaDiff) < 0.00864 * fabs(sinh(eta)) * etaSlice_)
0115                 continue;
0116               //if (sqrt(etaDiff*etaDiff + phiDiff*phiDiff) < 0.00864*fabs(sinh(eta))*intRadius_)
0117               //    continue;
0118               if ((etaDiff * etaDiff + phiDiff * phiDiff) < (0.000037325 * (cosh(2 * eta) - 1) * r2))
0119                 continue;
0120             }
0121           } else {
0122             if (fabs(etaDiff) < etaSlice_)
0123               continue;  // jurassic strip cut
0124             if (etaDiff * etaDiff + phiDiff * phiDiff < r2)
0125               continue;  // jurassic exclusion cone cut
0126           }
0127           //Check if RecHit is in SC
0128           if (vetoClustered_) {
0129             //Loop over basic clusters:
0130             bool isClustered = false;
0131             for (reco::CaloCluster_iterator bcIt = sc->clustersBegin(); bcIt != sc->clustersEnd(); ++bcIt) {
0132               for (rhIt = (*bcIt)->hitsAndFractions().begin(); rhIt != (*bcIt)->hitsAndFractions().end(); ++rhIt) {
0133                 if (rhIt->first == *i)
0134                   isClustered = true;
0135                 if (isClustered)
0136                   break;
0137               }
0138
0139               if (isClustered)
0140                 break;
0141             }  //end loop over basic clusters
0142
0143             if (isClustered)
0144               continue;
0145           }  //end if removeClustered
0146
0147           //std::cout << "detid " << j->detid() << std::endl;
0148           int severityFlag = ecalBarHits_ == nullptr ? -1 : sevLevel_->severityLevel(j->detid(), *ecalBarHits_);
0149           std::vector<int>::const_iterator sit =
0150               std::find(severitiesexcl_.begin(), severitiesexcl_.end(), severityFlag);
0151
0152           if (sit != severitiesexcl_.end())
0153             continue;
0154
0155           // new rechit flag checks
0156           //std::vector<int>::const_iterator vit = std::find(flags_.begin(),
0157           //                           flags_.end(),
0158           //                           j->recoFlag());
0159           //if (vit != flags_.end())
0160           //  continue;
0161           if (!j->checkFlag(EcalRecHit::kGood)) {
0162             if (j->checkFlags(flags_)) {
0163               continue;
0164             }
0165           }
0166
0167           float et = energy * std::sqrt(cell->getPosition().perp2() / cell->getPosition().mag2());
0168           if (et > etLow_ && energy > eLow_) {  //Changed energy --> fabs(energy) - now changed back to energy
0169             if (returnEt)
0170               energySum += et;
0171             else
0172               energySum += energy;
0173           }
0174
0175         }  //End if not end of list
0176       }  //End loop over rechits
0177     }  //End loop over barrel/endcap
0178   }  //End if caloHits_
0179
0180   return energySum;
0181 }
0182
0183 double EgammaRecHitIsolation::getSum_(const reco::SuperCluster* sc,
0184                                       bool returnEt,
0185                                       const EcalPFRecHitThresholds* thresholds) const {
0186   double energySum = 0.;
0187   if (!caloHits_.empty()) {
0188     //Take the SC position
0189
0190     const math::XYZPoint& theCaloPosition = sc->position();
0191     GlobalPoint pclu(theCaloPosition.x(), theCaloPosition.y(), theCaloPosition.z());
0192     double etaclus = pclu.eta();
0193     double phiclus = pclu.phi();
0194     double r2 = intRadius_ * intRadius_;
0195
0196     std::vector<std::pair<DetId, float> >::const_iterator rhIt;
0197
0198     for (int subdetnr = 0; subdetnr <= 1; subdetnr++) {  // look in barrel and endcap
0199       if (nullptr == subdet_[subdetnr])
0200         continue;
0201       CaloSubdetectorGeometry::DetIdSet chosen =
0202           subdet_[subdetnr]->getCells(pclu, extRadius_);  // select cells around cluster
0203       EcalRecHitCollection::const_iterator j = caloHits_.end();
0204       for (CaloSubdetectorGeometry::DetIdSet::const_iterator i = chosen.begin(); i != chosen.end();
0205            ++i) {  //loop selected cells
0206
0207         j = caloHits_.find(*i);      // find selected cell among rechits
0208         if (j != caloHits_.end()) {  // add rechit only if available
0209           const GlobalPoint& position = (theCaloGeom_.product())->getPosition(*i);
0210           double eta = position.eta();
0211           double phi = position.phi();
0212           double etaDiff = eta - etaclus;
0213           double phiDiff = reco::deltaPhi(phi, phiclus);
0214           double energy = j->energy();
0215
0216           float rhThres = (thresholds != nullptr) ? (*thresholds)[j->detid()] : 0.f;
0217           if (energy <= rhThres)
0218             continue;
0219
0220           if (useNumCrystals_) {
0221             if (fabs(etaclus) < 1.479) {  // Barrel num crystals, crystal width = 0.0174
0222               if (fabs(etaDiff) < 0.0174 * etaSlice_)
0223                 continue;
0224               // if ( sqrt(etaDiff*etaDiff + phiDiff*phiDiff) < 0.0174*intRadius_) continue;
0225               if ((etaDiff * etaDiff + phiDiff * phiDiff) < 0.00030276 * r2)
0226                 continue;
0227             } else {  // Endcap num crystals, crystal width = 0.00864*fabs(sinh(eta))
0228               if (fabs(etaDiff) < 0.00864 * fabs(sinh(eta)) * etaSlice_)
0229                 continue;
0230               // if ( sqrt(etaDiff*etaDiff + phiDiff*phiDiff) < 0.00864*fabs(sinh(eta))*intRadius_) continue;
0231               if ((etaDiff * etaDiff + phiDiff * phiDiff) < (0.000037325 * (cosh(2 * eta) - 1) * r2))
0232                 continue;
0233             }
0234           } else {
0235             if (fabs(etaDiff) < etaSlice_)
0236               continue;  // jurassic strip cut
0237             if (etaDiff * etaDiff + phiDiff * phiDiff < r2)
0238               continue;  // jurassic exclusion cone cut
0239           }
0240
0241           //Check if RecHit is in SC
0242           if (vetoClustered_) {
0243             //Loop over basic clusters:
0244             bool isClustered = false;
0245             for (reco::CaloCluster_iterator bcIt = sc->clustersBegin(); bcIt != sc->clustersEnd(); ++bcIt) {
0246               for (rhIt = (*bcIt)->hitsAndFractions().begin(); rhIt != (*bcIt)->hitsAndFractions().end(); ++rhIt) {
0247                 if (rhIt->first == *i)
0248                   isClustered = true;
0249                 if (isClustered)
0250                   break;
0251               }
0252               if (isClustered)
0253                 break;
0254             }  //end loop over basic clusters
0255
0256             if (isClustered)
0257               continue;
0258           }  //end if removeClustered
0259
0260           int severityFlag = sevLevel_->severityLevel(j->detid(), *ecalBarHits_);
0261           std::vector<int>::const_iterator sit =
0262               std::find(severitiesexcl_.begin(), severitiesexcl_.end(), severityFlag);
0263
0264           if (sit != severitiesexcl_.end())
0265             continue;
0266
0267           // new rechit flag checks
0268           //std::vector<int>::const_iterator vit = std::find(flags_.begin(),
0269           //                           flags_.end(),
0270           //                           j->recoFlag());
0271           //if (vit != flags_.end())
0272           //  continue;
0273           if (!j->checkFlag(EcalRecHit::kGood)) {
0274             if (j->checkFlags(flags_)) {
0275               continue;
0276             }
0277           }
0278
0279           double et = energy * position.perp() / position.mag();
0280           if (et > etLow_ && energy > eLow_) {  //Changed energy --> fabs(energy) -- then changed into energy
0281             if (returnEt)
0282               energySum += et;
0283             else
0284               energySum += energy;
0285           }
0286
0287         }  //End if not end of list
0288       }  //End loop over rechits
0289     }  //End loop over barrel/endcap
0290   }  //End if caloHits_
0291
0292   return energySum;
0293 }