File indexing completed on 2023-10-25 09:59:14
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, bool returnEt) const {
0068 double energySum = 0.;
0069 if (!caloHits_.empty()) {
0070
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++) {
0081 if (nullptr == subdet_[subdetnr])
0082 continue;
0083
0084 CaloSubdetectorGeometry::DetIdSet chosen =
0085 subdet_[subdetnr]->getCells(pclu, extRadius_);
0086 EcalRecHitCollection::const_iterator j = caloHits_.end();
0087
0088 for (CaloSubdetectorGeometry::DetIdSet::const_iterator i = chosen.begin(); i != chosen.end();
0089 ++i) {
0090 j = caloHits_.find(*i);
0091 if (j != caloHits_.end()) {
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) {
0101 if (fabs(etaDiff) < 0.0174 * etaSlice_)
0102 continue;
0103
0104
0105 if ((etaDiff * etaDiff + phiDiff * phiDiff) < 0.00030276 * r2)
0106 continue;
0107 } else {
0108 if (fabs(etaDiff) < 0.00864 * fabs(sinh(eta)) * etaSlice_)
0109 continue;
0110
0111
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;
0118 if (etaDiff * etaDiff + phiDiff * phiDiff < r2)
0119 continue;
0120 }
0121
0122 if (vetoClustered_) {
0123
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 }
0136
0137 if (isClustered)
0138 continue;
0139 }
0140
0141
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
0150
0151
0152
0153
0154
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_) {
0163 if (returnEt)
0164 energySum += et;
0165 else
0166 energySum += energy;
0167 }
0168
0169 }
0170 }
0171 }
0172 }
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
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++) {
0191 if (nullptr == subdet_[subdetnr])
0192 continue;
0193 CaloSubdetectorGeometry::DetIdSet chosen =
0194 subdet_[subdetnr]->getCells(pclu, extRadius_);
0195 EcalRecHitCollection::const_iterator j = caloHits_.end();
0196 for (CaloSubdetectorGeometry::DetIdSet::const_iterator i = chosen.begin(); i != chosen.end();
0197 ++i) {
0198
0199 j = caloHits_.find(*i);
0200 if (j != caloHits_.end()) {
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) {
0210 if (fabs(etaDiff) < 0.0174 * etaSlice_)
0211 continue;
0212
0213 if ((etaDiff * etaDiff + phiDiff * phiDiff) < 0.00030276 * r2)
0214 continue;
0215 } else {
0216 if (fabs(etaDiff) < 0.00864 * fabs(sinh(eta)) * etaSlice_)
0217 continue;
0218
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;
0225 if (etaDiff * etaDiff + phiDiff * phiDiff < r2)
0226 continue;
0227 }
0228
0229
0230 if (vetoClustered_) {
0231
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 }
0243
0244 if (isClustered)
0245 continue;
0246 }
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
0256
0257
0258
0259
0260
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_) {
0269 if (returnEt)
0270 energySum += et;
0271 else
0272 energySum += energy;
0273 }
0274
0275 }
0276 }
0277 }
0278 }
0279
0280 return energySum;
0281 }