File indexing completed on 2023-03-17 11:17:43
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 &)) const {
0130 const HcalDetId hid(hit.detid());
0131 const int hd = hid.depth(), he = hid.ieta(), hp = hid.iphi();
0132 const int h1 = hd - 1;
0133
0134 if (include_or_exclude == -1 and (he != ieta or hp != iphi))
0135 return 0.;
0136
0137 if (include_or_exclude == 1 and (he == ieta and hp == iphi))
0138 return 0.;
0139
0140 if ((hid.subdet() == HcalBarrel and (hd < 1 or hd > int(eThresHB_.size()))) or
0141 (hid.subdet() == HcalEndcap and (hd < 1 or hd > int(eThresHE_.size()))))
0142 edm::LogWarning("EgammaHcalIsolation")
0143 << " hit in subdet " << hid.subdet() << " has an unaccounted for depth of " << hd << "!!";
0144
0145 const bool right_depth = (depth == 0 or hd == depth);
0146 if (!right_depth)
0147 return 0.;
0148
0149 const bool goodHBe = hid.subdet() == HcalBarrel and hit.energy() > eThresHB_[h1];
0150 const bool goodHEe = hid.subdet() == HcalEndcap and hit.energy() > eThresHE_[h1];
0151 if (!(goodHBe or goodHEe))
0152 return 0.;
0153
0154 const auto phit = caloGeometry_.getGeometry(hit.detid())->repPos();
0155 const float phitEta = phit.eta();
0156
0157 if (extIncRule_ == InclusionRule::withinConeAroundCluster or intIncRule_ == InclusionRule::withinConeAroundCluster) {
0158 auto const dR2 = deltaR2(pcluEta, pcluPhi, phitEta, phit.phi());
0159 if ((extIncRule_ == InclusionRule::withinConeAroundCluster and dR2 > extRadius_) or
0160 (intIncRule_ == InclusionRule::withinConeAroundCluster and dR2 < intRadius_))
0161 return 0.;
0162 }
0163
0164 DetId did = hcalTopology_.idFront(hid);
0165 const uint32_t flag = hit.flags();
0166 const uint32_t dbflag = hcalChStatus_.getValues(did)->getValue();
0167 int severity = hcalSevLvlComputer_.getSeverityLevel(did, flag, dbflag);
0168 bool recovered = hcalSevLvlComputer_.recoveredRecHit(did, flag);
0169
0170 const double het = hit.energy() * scaleToEt(phitEta);
0171 const bool goodHB = goodHBe and (severity <= maxSeverityHB_ or recovered) and het > etThresHB_[h1];
0172 const bool goodHE = goodHEe and (severity <= maxSeverityHE_ or recovered) and het > etThresHE_[h1];
0173
0174 if (goodHB or goodHE)
0175 return hit.energy() * scale(phitEta);
0176
0177 return 0.;
0178 }
0179
0180 double EgammaHcalIsolation::getHcalSum(const GlobalPoint &pclu,
0181 int depth,
0182 int ieta,
0183 int iphi,
0184 int include_or_exclude,
0185 double (*scale)(const double &)) const {
0186 double sum = 0.;
0187 const float pcluEta = pclu.eta();
0188 const float pcluPhi = pclu.phi();
0189 for (const auto &hit : mhbhe_)
0190 sum += goodHitEnergy(pcluEta, pcluPhi, hit, depth, ieta, iphi, include_or_exclude, scale);
0191
0192 return sum;
0193 }
