File indexing completed on 2024-04-06 12:24:56
0001
0002
0003
0004
0005
0006
0007
0008
0009 #include <Math/VectorUtil.h>
0010
0011
0012 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0013 #include "DataFormats/HcalDetId/interface/HcalDetId.h"
0014 #include "DataFormats/Math/interface/deltaR.h"
0015 #include "DataFormats/RecoCandidate/interface/RecoCandidate.h"
0016 #include "Geometry/CaloGeometry/interface/CaloCellGeometry.h"
0017 #include "Geometry/CommonDetUnit/interface/TrackingGeometry.h"
0018 #include "Geometry/Records/interface/TrackerDigiGeometryRecord.h"
0019 #include "RecoEgamma/EgammaIsolationAlgos/interface/EgammaHcalIsolation.h"
0020
0021 double scaleToE(const double &eta) { return 1.; }
0022 double scaleToEt(const double &eta) { return std::sin(2. * std::atan(std::exp(-eta))); }
0023
0024 EgammaHcalIsolation::EgammaHcalIsolation(InclusionRule extIncRule,
0025 double extRadius,
0026 InclusionRule intIncRule,
0027 double intRadius,
0028 const arrayHB &eThresHB,
0029 const arrayHB &etThresHB,
0030 int maxSeverityHB,
0031 const arrayHE &eThresHE,
0032 const arrayHE &etThresHE,
0033 int maxSeverityHE,
0034 const HBHERecHitCollection &mhbhe,
0035 edm::ESHandle<CaloGeometry> caloGeometry,
0036 edm::ESHandle<HcalTopology> hcalTopology,
0037 edm::ESHandle<HcalChannelQuality> hcalChStatus,
0038 edm::ESHandle<HcalSeverityLevelComputer> hcalSevLvlComputer,
0039 edm::ESHandle<CaloTowerConstituentsMap> towerMap)
0040 : extIncRule_(extIncRule),
0041 extRadius_(extRadius * extRadius),
0042 intIncRule_(intIncRule),
0043 intRadius_(intRadius * intRadius),
0044 maxSeverityHB_(maxSeverityHB),
0045 maxSeverityHE_(maxSeverityHE),
0046 mhbhe_(mhbhe),
0047 caloGeometry_(*caloGeometry.product()),
0048 hcalTopology_(*hcalTopology.product()),
0049 hcalChStatus_(*hcalChStatus.product()),
0050 hcalSevLvlComputer_(*hcalSevLvlComputer.product()),
0051 towerMap_(*towerMap.product()) {
0052 eThresHB_ = eThresHB;
0053 etThresHB_ = etThresHB;
0054 eThresHE_ = eThresHE;
0055 etThresHE_ = etThresHE;
0056
0057
0058 if (extIncRule_ == InclusionRule::isBehindClusterSeed and intIncRule_ == InclusionRule::withinConeAroundCluster) {
0059 extRadius_ = 0.;
0060 intRadius_ = 0.;
0061 } else if (extIncRule_ == InclusionRule::withinConeAroundCluster and
0062 intIncRule_ == InclusionRule::isBehindClusterSeed) {
0063 intRadius_ = 0.;
0064 } else if (extIncRule_ == InclusionRule::isBehindClusterSeed and intIncRule_ == InclusionRule::isBehindClusterSeed) {
0065 edm::LogWarning("EgammaHcalIsolation")
0066 << " external and internal rechit inclusion rules can't both be isBehindClusterSeed."
0067 << " Setting both to withinConeAroundCluster!";
0068 extIncRule_ = InclusionRule::withinConeAroundCluster;
0069 intIncRule_ = InclusionRule::withinConeAroundCluster;
0070 }
0071 }
0072
0073 EgammaHcalIsolation::EgammaHcalIsolation(InclusionRule extIncRule,
0074 double extRadius,
0075 InclusionRule intIncRule,
0076 double intRadius,
0077 const arrayHB &eThresHB,
0078 const arrayHB &etThresHB,
0079 int maxSeverityHB,
0080 const arrayHE &eThresHE,
0081 const arrayHE &etThresHE,
0082 int maxSeverityHE,
0083 const HBHERecHitCollection &mhbhe,
0084 const CaloGeometry &caloGeometry,
0085 const HcalTopology &hcalTopology,
0086 const HcalChannelQuality &hcalChStatus,
0087 const HcalSeverityLevelComputer &hcalSevLvlComputer,
0088 const CaloTowerConstituentsMap &towerMap)
0089 : extIncRule_(extIncRule),
0090 extRadius_(extRadius * extRadius),
0091 intIncRule_(intIncRule),
0092 intRadius_(intRadius * intRadius),
0093 maxSeverityHB_(maxSeverityHB),
0094 maxSeverityHE_(maxSeverityHE),
0095 mhbhe_(mhbhe),
0096 caloGeometry_(caloGeometry),
0097 hcalTopology_(hcalTopology),
0098 hcalChStatus_(hcalChStatus),
0099 hcalSevLvlComputer_(hcalSevLvlComputer),
0100 towerMap_(towerMap) {
0101 eThresHB_ = eThresHB;
0102 etThresHB_ = etThresHB;
0103 eThresHE_ = eThresHE;
0104 etThresHE_ = etThresHE;
0105
0106
0107 if (extIncRule_ == InclusionRule::isBehindClusterSeed and intIncRule_ == InclusionRule::withinConeAroundCluster) {
0108 extRadius_ = 0.;
0109 intRadius_ = 0.;
0110 } else if (extIncRule_ == InclusionRule::withinConeAroundCluster and
0111 intIncRule_ == InclusionRule::isBehindClusterSeed) {
0112 intRadius_ = 0.;
0113 } else if (extIncRule_ == InclusionRule::isBehindClusterSeed and intIncRule_ == InclusionRule::isBehindClusterSeed) {
0114 edm::LogWarning("EgammaHcalIsolation")
0115 << " external and internal rechit inclusion rules can't both be isBehindClusterSeed."
0116 << " Setting both to withinConeAroundCluster!";
0117 extIncRule_ = InclusionRule::withinConeAroundCluster;
0118 intIncRule_ = InclusionRule::withinConeAroundCluster;
0119 }
0120 }
0121
0122 double EgammaHcalIsolation::goodHitEnergy(float pcluEta,
0123 float pcluPhi,
0124 const HBHERecHit &hit,
0125 int depth,
0126 int ieta,
0127 int iphi,
0128 int include_or_exclude,
0129 double (*scale)(const double &),
0130 const HcalPFCuts *hcalCuts) const {
0131 const HcalDetId hid(hit.detid());
0132 const int hd = hid.depth(), he = hid.ieta(), hp = hid.iphi();
0133 const int h1 = hd - 1;
0134 double thresholdE = 0.;
0135
0136 if (include_or_exclude == -1 and (he != ieta or hp != iphi))
0137 return 0.;
0138
0139 if (include_or_exclude == 1 and (he == ieta and hp == iphi))
0140 return 0.;
0141
0142 if ((hid.subdet() == HcalBarrel and (hd < 1 or hd > int(eThresHB_.size()))) or
0143 (hid.subdet() == HcalEndcap and (hd < 1 or hd > int(eThresHE_.size()))))
0144 edm::LogWarning("EgammaHcalIsolation")
0145 << " hit in subdet " << hid.subdet() << " has an unaccounted for depth of " << hd << "!!";
0146
0147 const bool right_depth = (depth == 0 or hd == depth);
0148 if (!right_depth)
0149 return 0.;
0150
0151 bool goodHBe = hid.subdet() == HcalBarrel and hit.energy() > eThresHB_[h1];
0152 bool goodHEe = hid.subdet() == HcalEndcap and hit.energy() > eThresHE_[h1];
0153
0154 if (hcalCuts != nullptr) {
0155 const HcalPFCut *cutValue = hcalCuts->getValues(hid.rawId());
0156 thresholdE = cutValue->noiseThreshold();
0157 goodHBe = hid.subdet() == HcalBarrel and hit.energy() > thresholdE;
0158 goodHEe = hid.subdet() == HcalEndcap and hit.energy() > thresholdE;
0159 }
0160
0161 if (!(goodHBe or goodHEe))
0162 return 0.;
0163
0164 const auto phit = caloGeometry_.getGeometry(hit.detid())->repPos();
0165 const float phitEta = phit.eta();
0166
0167 if (extIncRule_ == InclusionRule::withinConeAroundCluster or intIncRule_ == InclusionRule::withinConeAroundCluster) {
0168 auto const dR2 = deltaR2(pcluEta, pcluPhi, phitEta, phit.phi());
0169 if ((extIncRule_ == InclusionRule::withinConeAroundCluster and dR2 > extRadius_) or
0170 (intIncRule_ == InclusionRule::withinConeAroundCluster and dR2 < intRadius_))
0171 return 0.;
0172 }
0173
0174 DetId did = hcalTopology_.idFront(hid);
0175 const uint32_t flag = hit.flags();
0176 const uint32_t dbflag = hcalChStatus_.getValues(did)->getValue();
0177 int severity = hcalSevLvlComputer_.getSeverityLevel(did, flag, dbflag);
0178 bool recovered = hcalSevLvlComputer_.recoveredRecHit(did, flag);
0179
0180 const double het = hit.energy() * scaleToEt(phitEta);
0181 const bool goodHB = goodHBe and (severity <= maxSeverityHB_ or recovered) and het > etThresHB_[h1];
0182 const bool goodHE = goodHEe and (severity <= maxSeverityHE_ or recovered) and het > etThresHE_[h1];
0183
0184 if (goodHB or goodHE)
0185 return hit.energy() * scale(phitEta);
0186
0187 return 0.;
0188 }
0189
0190 double EgammaHcalIsolation::getHcalSum(const GlobalPoint &pclu,
0191 int depth,
0192 int ieta,
0193 int iphi,
0194 int include_or_exclude,
0195 double (*scale)(const double &),
0196 const HcalPFCuts *hcalCuts) const {
0197 double sum = 0.;
0198 const float pcluEta = pclu.eta();
0199 const float pcluPhi = pclu.phi();
0200 for (const auto &hit : mhbhe_)
0201 sum += goodHitEnergy(pcluEta, pcluPhi, hit, depth, ieta, iphi, include_or_exclude, scale, hcalCuts);
0202
0203 return sum;
0204 }
