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, bool returnEt) const {
0068   double energySum = 0.;
0069   if (!caloHits_.empty()) {
0070     //Take the SC position
0071     reco::SuperClusterRef sc = emObject->get<reco::SuperClusterRef>();
0072     math::XYZPoint const& theCaloPosition = sc.get()->position();
0073     GlobalPoint pclu(theCaloPosition.x(), theCaloPosition.y(), theCaloPosition.z());
0074     float etaclus = pclu.eta();
0075     float phiclus = pclu.phi();
0076     float r2 = intRadius_ * intRadius_;
0077
0078     std::vector<std::pair<DetId, float> >::const_iterator rhIt;
0079
0080     for (int subdetnr = 0; subdetnr <= 1; subdetnr++) {  // look in barrel and endcap
0081       if (nullptr == subdet_[subdetnr])
0082         continue;
0083
0084       CaloSubdetectorGeometry::DetIdSet chosen =
0085           subdet_[subdetnr]->getCells(pclu, extRadius_);  // select cells around cluster
0086       EcalRecHitCollection::const_iterator j = caloHits_.end();
0087
0088       for (CaloSubdetectorGeometry::DetIdSet::const_iterator i = chosen.begin(); i != chosen.end();
0089            ++i) {                    //loop selected cells
0090         j = caloHits_.find(*i);      // find selected cell among rechits
0091         if (j != caloHits_.end()) {  // add rechit only if available
0092           auto cell = theCaloGeom_->getGeometry(*i);
0093           float eta = cell->etaPos();
0094           float phi = cell->phiPos();
0095           float etaDiff = eta - etaclus;
0096           float phiDiff = reco::deltaPhi(phi, phiclus);
0097           float energy = j->energy();
0098
0099           if (useNumCrystals_) {
0100             if (fabs(etaclus) < 1.479) {  // Barrel num crystals, crystal width = 0.0174
0101               if (fabs(etaDiff) < 0.0174 * etaSlice_)
0102                 continue;
0103               //if (sqrt(etaDiff*etaDiff + phiDiff*phiDiff) < 0.0174*intRadius_)
0104               //continue;
0105               if ((etaDiff * etaDiff + phiDiff * phiDiff) < 0.00030276 * r2)
0106                 continue;
0107             } else {  // Endcap num crystals, crystal width = 0.00864*fabs(sinh(eta))
0108               if (fabs(etaDiff) < 0.00864 * fabs(sinh(eta)) * etaSlice_)
0109                 continue;
0110               //if (sqrt(etaDiff*etaDiff + phiDiff*phiDiff) < 0.00864*fabs(sinh(eta))*intRadius_)
0111               //    continue;
0112               if ((etaDiff * etaDiff + phiDiff * phiDiff) < (0.000037325 * (cosh(2 * eta) - 1) * r2))
0113                 continue;
0114             }
0115           } else {
0116             if (fabs(etaDiff) < etaSlice_)
0117               continue;  // jurassic strip cut
0118             if (etaDiff * etaDiff + phiDiff * phiDiff < r2)
0119               continue;  // jurassic exclusion cone cut
0120           }
0121           //Check if RecHit is in SC
0122           if (vetoClustered_) {
0123             //Loop over basic clusters:
0124             bool isClustered = false;
0125             for (reco::CaloCluster_iterator bcIt = sc->clustersBegin(); bcIt != sc->clustersEnd(); ++bcIt) {
0126               for (rhIt = (*bcIt)->hitsAndFractions().begin(); rhIt != (*bcIt)->hitsAndFractions().end(); ++rhIt) {
0127                 if (rhIt->first == *i)
0128                   isClustered = true;
0129                 if (isClustered)
0130                   break;
0131               }
0132
0133               if (isClustered)
0134                 break;
0135             }  //end loop over basic clusters
0136
0137             if (isClustered)
0138               continue;
0139           }  //end if removeClustered
0140
0141           //std::cout << "detid " << j->detid() << std::endl;
0142           int severityFlag = ecalBarHits_ == nullptr ? -1 : sevLevel_->severityLevel(j->detid(), *ecalBarHits_);
0143           std::vector<int>::const_iterator sit =
0144               std::find(severitiesexcl_.begin(), severitiesexcl_.end(), severityFlag);
0145
0146           if (sit != severitiesexcl_.end())
0147             continue;
0148
0149           // new rechit flag checks
0150           //std::vector<int>::const_iterator vit = std::find(flags_.begin(),
0151           //                           flags_.end(),
0152           //                           j->recoFlag());
0153           //if (vit != flags_.end())
0154           //  continue;
0155           if (!j->checkFlag(EcalRecHit::kGood)) {
0156             if (j->checkFlags(flags_)) {
0157               continue;
0158             }
0159           }
0160
0161           float et = energy * std::sqrt(cell->getPosition().perp2() / cell->getPosition().mag2());
0162           if (et > etLow_ && energy > eLow_) {  //Changed energy --> fabs(energy) - now changed back to energy
0163             if (returnEt)
0164               energySum += et;
0165             else
0166               energySum += energy;
0167           }
0168
0169         }  //End if not end of list
0170       }    //End loop over rechits
0171     }      //End loop over barrel/endcap
0172   }        //End if caloHits_
0173
0174   return energySum;
0175 }
0176
0177 double EgammaRecHitIsolation::getSum_(const reco::SuperCluster* sc, bool returnEt) const {
0178   double energySum = 0.;
0179   if (!caloHits_.empty()) {
0180     //Take the SC position
0181
0182     const math::XYZPoint& theCaloPosition = sc->position();
0183     GlobalPoint pclu(theCaloPosition.x(), theCaloPosition.y(), theCaloPosition.z());
0184     double etaclus = pclu.eta();
0185     double phiclus = pclu.phi();
0186     double r2 = intRadius_ * intRadius_;
0187
0188     std::vector<std::pair<DetId, float> >::const_iterator rhIt;
0189
0190     for (int subdetnr = 0; subdetnr <= 1; subdetnr++) {  // look in barrel and endcap
0191       if (nullptr == subdet_[subdetnr])
0192         continue;
0193       CaloSubdetectorGeometry::DetIdSet chosen =
0194           subdet_[subdetnr]->getCells(pclu, extRadius_);  // select cells around cluster
0195       EcalRecHitCollection::const_iterator j = caloHits_.end();
0196       for (CaloSubdetectorGeometry::DetIdSet::const_iterator i = chosen.begin(); i != chosen.end();
0197            ++i) {  //loop selected cells
0198
0199         j = caloHits_.find(*i);      // find selected cell among rechits
0200         if (j != caloHits_.end()) {  // add rechit only if available
0201           const GlobalPoint& position = (theCaloGeom_.product())->getPosition(*i);
0202           double eta = position.eta();
0203           double phi = position.phi();
0204           double etaDiff = eta - etaclus;
0205           double phiDiff = reco::deltaPhi(phi, phiclus);
0206           double energy = j->energy();
0207
0208           if (useNumCrystals_) {
0209             if (fabs(etaclus) < 1.479) {  // Barrel num crystals, crystal width = 0.0174
0210               if (fabs(etaDiff) < 0.0174 * etaSlice_)
0211                 continue;
0212               // if ( sqrt(etaDiff*etaDiff + phiDiff*phiDiff) < 0.0174*intRadius_) continue;
0213               if ((etaDiff * etaDiff + phiDiff * phiDiff) < 0.00030276 * r2)
0214                 continue;
0215             } else {  // Endcap num crystals, crystal width = 0.00864*fabs(sinh(eta))
0216               if (fabs(etaDiff) < 0.00864 * fabs(sinh(eta)) * etaSlice_)
0217                 continue;
0218               // if ( sqrt(etaDiff*etaDiff + phiDiff*phiDiff) < 0.00864*fabs(sinh(eta))*intRadius_) continue;
0219               if ((etaDiff * etaDiff + phiDiff * phiDiff) < (0.000037325 * (cosh(2 * eta) - 1) * r2))
0220                 continue;
0221             }
0222           } else {
0223             if (fabs(etaDiff) < etaSlice_)
0224               continue;  // jurassic strip cut
0225             if (etaDiff * etaDiff + phiDiff * phiDiff < r2)
0226               continue;  // jurassic exclusion cone cut
0227           }
0228
0229           //Check if RecHit is in SC
0230           if (vetoClustered_) {
0231             //Loop over basic clusters:
0232             bool isClustered = false;
0233             for (reco::CaloCluster_iterator bcIt = sc->clustersBegin(); bcIt != sc->clustersEnd(); ++bcIt) {
0234               for (rhIt = (*bcIt)->hitsAndFractions().begin(); rhIt != (*bcIt)->hitsAndFractions().end(); ++rhIt) {
0235                 if (rhIt->first == *i)
0236                   isClustered = true;
0237                 if (isClustered)
0238                   break;
0239               }
0240               if (isClustered)
0241                 break;
0242             }  //end loop over basic clusters
0243
0244             if (isClustered)
0245               continue;
0246           }  //end if removeClustered
0247
0248           int severityFlag = sevLevel_->severityLevel(j->detid(), *ecalBarHits_);
0249           std::vector<int>::const_iterator sit =
0250               std::find(severitiesexcl_.begin(), severitiesexcl_.end(), severityFlag);
0251
0252           if (sit != severitiesexcl_.end())
0253             continue;
0254
0255           // new rechit flag checks
0256           //std::vector<int>::const_iterator vit = std::find(flags_.begin(),
0257           //                           flags_.end(),
0258           //                           j->recoFlag());
0259           //if (vit != flags_.end())
0260           //  continue;
0261           if (!j->checkFlag(EcalRecHit::kGood)) {
0262             if (j->checkFlags(flags_)) {
0263               continue;
0264             }
0265           }
0266
0267           double et = energy * position.perp() / position.mag();
0268           if (et > etLow_ && energy > eLow_) {  //Changed energy --> fabs(energy) -- then changed into energy
0269             if (returnEt)
0270               energySum += et;
0271             else
0272               energySum += energy;
0273           }
0274
0275         }  //End if not end of list
0276       }    //End loop over rechits
0277     }      //End loop over barrel/endcap
0278   }        //End if caloHits_
0279
0280   return energySum;
0281 }