File indexing completed on 2024-09-07 04:37:28
0001
0002
0003
0004
0005
0006
0007
0008
0009 #include <vector>
0010 #include <functional>
0011
0012
0013 #include <Math/VectorUtil.h>
0014
0015
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 :
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
0054 severitiesexcl_(0),
0055
0056
0057
0058 flags_(0) {
0059
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
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++) {
0083 if (nullptr == subdet_[subdetnr])
0084 continue;
0085
0086 CaloSubdetectorGeometry::DetIdSet chosen =
0087 subdet_[subdetnr]->getCells(pclu, extRadius_);
0088 EcalRecHitCollection::const_iterator j = caloHits_.end();
0089
0090 for (CaloSubdetectorGeometry::DetIdSet::const_iterator i = chosen.begin(); i != chosen.end();
0091 ++i) {
0092 j = caloHits_.find(*i);
0093 if (j != caloHits_.end()) {
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) {
0107 if (fabs(etaDiff) < 0.0174 * etaSlice_)
0108 continue;
0109
0110
0111 if ((etaDiff * etaDiff + phiDiff * phiDiff) < 0.00030276 * r2)
0112 continue;
0113 } else {
0114 if (fabs(etaDiff) < 0.00864 * fabs(sinh(eta)) * etaSlice_)
0115 continue;
0116
0117
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;
0124 if (etaDiff * etaDiff + phiDiff * phiDiff < r2)
0125 continue;
0126 }
0127
0128 if (vetoClustered_) {
0129
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 }
0142
0143 if (isClustered)
0144 continue;
0145 }
0146
0147
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
0156
0157
0158
0159
0160
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_) {
0169 if (returnEt)
0170 energySum += et;
0171 else
0172 energySum += energy;
0173 }
0174
0175 }
0176 }
0177 }
0178 }
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
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++) {
0199 if (nullptr == subdet_[subdetnr])
0200 continue;
0201 CaloSubdetectorGeometry::DetIdSet chosen =
0202 subdet_[subdetnr]->getCells(pclu, extRadius_);
0203 EcalRecHitCollection::const_iterator j = caloHits_.end();
0204 for (CaloSubdetectorGeometry::DetIdSet::const_iterator i = chosen.begin(); i != chosen.end();
0205 ++i) {
0206
0207 j = caloHits_.find(*i);
0208 if (j != caloHits_.end()) {
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) {
0222 if (fabs(etaDiff) < 0.0174 * etaSlice_)
0223 continue;
0224
0225 if ((etaDiff * etaDiff + phiDiff * phiDiff) < 0.00030276 * r2)
0226 continue;
0227 } else {
0228 if (fabs(etaDiff) < 0.00864 * fabs(sinh(eta)) * etaSlice_)
0229 continue;
0230
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;
0237 if (etaDiff * etaDiff + phiDiff * phiDiff < r2)
0238 continue;
0239 }
0240
0241
0242 if (vetoClustered_) {
0243
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 }
0255
0256 if (isClustered)
0257 continue;
0258 }
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
0268
0269
0270
0271
0272
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_) {
0281 if (returnEt)
0282 energySum += et;
0283 else
0284 energySum += energy;
0285 }
0286
0287 }
0288 }
0289 }
0290 }
0291
0292 return energySum;
0293 }