File indexing completed on 2024-04-06 11:59:20
0001 #include "Calibration/IsolatedParticles/interface/eECALMatrix.h"
0002 #include "Calibration/IsolatedParticles/interface/FindCaloHit.h"
0003 #include "DataFormats/CaloTowers/interface/CaloTowerDetId.h"
0004 #include "DataFormats/EcalDetId/interface/EBDetId.h"
0005 #include "DataFormats/EcalDetId/interface/EEDetId.h"
0006 #include "DataFormats/EcalDetId/interface/ESDetId.h"
0007 #include "DataFormats/EcalDetId/interface/EcalSubdetector.h"
0008 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0009 #include "RecoLocalCalo/EcalRecAlgos/interface/EcalSeverityLevelAlgo.h"
0010
0011 #include <sstream>
0012
0013 namespace spr {
0014
0015 std::pair<double, bool> energyECAL(const DetId& id,
0016 edm::Handle<EcalRecHitCollection>& hitsEC,
0017 const EcalSeverityLevelAlgo* sevlv,
0018 bool testSpike,
0019 double tMin,
0020 double tMax,
0021 bool debug) {
0022 std::vector<EcalRecHitCollection::const_iterator> hits;
0023 spr::findHit(hitsEC, id, hits, debug);
0024
0025 std::ostringstream st1;
0026 if (debug)
0027 st1 << "Xtal 0x" << std::hex << id() << std::dec;
0028 const EcalRecHitCollection* recHitsEC = (hitsEC.isValid()) ? hitsEC.product() : nullptr;
0029 bool flag = (!testSpike) ? true : (sevlv->severityLevel(id, (*recHitsEC)) != EcalSeverityLevel::kWeird);
0030 double ener(0);
0031 for (const auto& hit : hits) {
0032 double en(0), tt(0);
0033 if (hit != hitsEC->end()) {
0034 en = hit->energy();
0035 tt = hit->time();
0036 }
0037 if (debug)
0038 st1 << " " << tt << " " << en;
0039 if (tt > tMin && tt < tMax)
0040 ener += en;
0041 }
0042 if (debug) {
0043 if (!flag)
0044 st1 << " detected to be a spike";
0045 edm::LogVerbatim("IsoTrack") << st1.str();
0046 }
0047 return std::pair<double, bool>(ener, flag);
0048 }
0049
0050 std::pair<double, bool> energyECAL(const std::vector<DetId>& vdets,
0051 edm::Handle<EcalRecHitCollection>& hitsEC,
0052 const EcalSeverityLevelAlgo* sevlv,
0053 bool noThrCut,
0054 bool testSpike,
0055 double eThr,
0056 double tMin,
0057 double tMax,
0058 bool debug) {
0059 bool flag(true);
0060 double energySum(0.0);
0061 for (const auto& id : vdets) {
0062 if (id != DetId(0)) {
0063 std::pair<double, bool> ecalEn = spr::energyECAL(id, hitsEC, sevlv, testSpike, tMin, tMax, debug);
0064 if (!ecalEn.second)
0065 flag = false;
0066 if ((ecalEn.first > eThr) || noThrCut)
0067 energySum += ecalEn.first;
0068 }
0069 }
0070 if (debug)
0071 edm::LogVerbatim("IsoTrack") << "energyECAL: energySum = " << energySum << " flag = " << flag;
0072 return std::pair<double, bool>(energySum, flag);
0073 }
0074
0075 std::pair<double, bool> eECALmatrix(const DetId& detId,
0076 edm::Handle<EcalRecHitCollection>& hitsEB,
0077 edm::Handle<EcalRecHitCollection>& hitsEE,
0078 const EcalChannelStatus& chStatus,
0079 const CaloGeometry* geo,
0080 const CaloTopology* caloTopology,
0081 const EcalSeverityLevelAlgo* sevlv,
0082 int ieta,
0083 int iphi,
0084 double ebThr,
0085 double eeThr,
0086 double tMin,
0087 double tMax,
0088 bool debug) {
0089 std::vector<DetId> vdets;
0090 spr::matrixECALIds(detId, ieta, iphi, geo, caloTopology, vdets, debug);
0091 if (debug)
0092 edm::LogVerbatim("IsoTrack") << "Inside eECALmatrix " << 2 * ieta + 1 << "X" << 2 * iphi + 1 << " nXtals "
0093 << vdets.size();
0094 bool flag(true);
0095 for (const auto& id : vdets) {
0096 if ((id.det() == DetId::Ecal) && (id.subdetId() == EcalBarrel)) {
0097 if (sevlv->severityLevel(id, (*hitsEB)) == EcalSeverityLevel::kWeird)
0098 flag = false;
0099 } else if ((id.det() == DetId::Ecal) && (id.subdetId() == EcalEndcap)) {
0100 if (sevlv->severityLevel(id, (*hitsEE)) == EcalSeverityLevel::kWeird)
0101 flag = false;
0102 }
0103 }
0104 return std::pair<double, bool>(spr::energyECAL(vdets, hitsEB, hitsEE, ebThr, eeThr, tMin, tMax, debug), flag);
0105 }
0106
0107 std::pair<double, bool> eECALmatrix(const DetId& detId,
0108 edm::Handle<EcalRecHitCollection>& hitsEB,
0109 edm::Handle<EcalRecHitCollection>& hitsEE,
0110 const EcalChannelStatus& chStatus,
0111 const CaloGeometry* geo,
0112 const CaloTopology* caloTopology,
0113 const EcalSeverityLevelAlgo* sevlv,
0114 const EcalTrigTowerConstituentsMap& ttMap,
0115 int ieta,
0116 int iphi,
0117 double ebThr,
0118 double eeThr,
0119 double tMin,
0120 double tMax,
0121 bool debug) {
0122 std::vector<DetId> vdets;
0123 spr::matrixECALIds(detId, ieta, iphi, geo, caloTopology, vdets, debug);
0124 if (debug)
0125 edm::LogVerbatim("IsoTrack") << "Inside eECALmatrix " << 2 * ieta + 1 << "X" << 2 * iphi + 1 << " nXtals "
0126 << vdets.size();
0127
0128 bool flag(true);
0129 double energySum = 0.0;
0130 for (const auto& id : vdets) {
0131 if ((id != DetId(0)) && (id.det() == DetId::Ecal) &&
0132 ((id.subdetId() == EcalBarrel) || (id.subdetId() == EcalEndcap))) {
0133 double eTower = spr::energyECALTower(id, hitsEB, hitsEE, ttMap, debug);
0134 bool ok = (id.subdetId() == EcalBarrel) ? (eTower > ebThr) : (eTower > eeThr);
0135 if (debug && (!ok))
0136 edm::LogVerbatim("IsoTrack") << "Crystal 0x" << std::hex << id() << std::dec << " Flag " << ok;
0137 if (ok) {
0138 std::pair<double, bool> ecalEn = (id.subdetId() == EcalBarrel)
0139 ? spr::energyECAL(id, hitsEB, sevlv, true, tMin, tMax, debug)
0140 : spr::energyECAL(id, hitsEE, sevlv, false, tMin, tMax, debug);
0141 if (!ecalEn.second)
0142 flag = false;
0143 energySum += ecalEn.first;
0144 }
0145 }
0146 }
0147 if (debug)
0148 edm::LogVerbatim("IsoTrack") << "energyECAL: energySum = " << energySum << " flag = " << flag;
0149 return std::pair<double, bool>(energySum, flag);
0150 }
0151
0152 std::pair<double, bool> eECALmatrix(const HcalDetId& detId,
0153 edm::Handle<EcalRecHitCollection>& hitsEB,
0154 edm::Handle<EcalRecHitCollection>& hitsEE,
0155 const CaloGeometry* geo,
0156 const CaloTowerConstituentsMap* ctmap,
0157 const EcalSeverityLevelAlgo* sevlv,
0158 double ebThr,
0159 double eeThr,
0160 double tMin,
0161 double tMax,
0162 bool debug) {
0163 CaloTowerDetId tower = ctmap->towerOf(detId);
0164 std::vector<DetId> ids = ctmap->constituentsOf(tower);
0165 if (debug) {
0166 edm::LogVerbatim("IsoTrack") << "eECALmatrix: " << detId << " belongs to " << tower << " which has " << ids.size()
0167 << " constituents";
0168 for (unsigned int i = 0; i < ids.size(); ++i) {
0169 std::ostringstream st1;
0170 st1 << "[" << i << "] " << std::hex << ids[i].rawId() << std::dec;
0171 if (ids[i].det() == DetId::Ecal && ids[i].subdetId() == EcalBarrel) {
0172 st1 << " " << EBDetId(ids[i]);
0173 } else if (ids[i].det() == DetId::Ecal && ids[i].subdetId() == EcalEndcap) {
0174 st1 << " " << EEDetId(ids[i]);
0175 } else if (ids[i].det() == DetId::Ecal && ids[i].subdetId() == EcalPreshower) {
0176 st1 << " " << ESDetId(ids[i]);
0177 } else if (ids[i].det() == DetId::Hcal) {
0178 st1 << " " << HcalDetId(ids[i]);
0179 }
0180 edm::LogVerbatim("IsoTrack") << st1.str();
0181 }
0182 }
0183
0184 std::vector<DetId> idEBEE;
0185 bool flag(true);
0186 for (const auto& id : ids) {
0187 if ((id.det() == DetId::Ecal) && (id.subdetId() == EcalBarrel)) {
0188 idEBEE.emplace_back(id);
0189 if (sevlv->severityLevel(id, (*hitsEB)) == EcalSeverityLevel::kWeird)
0190 flag = false;
0191 } else if ((id.det() == DetId::Ecal) && (id.subdetId() == EcalEndcap)) {
0192 idEBEE.emplace_back(id);
0193 if (sevlv->severityLevel(id, (*hitsEE)) == EcalSeverityLevel::kWeird)
0194 flag = false;
0195 }
0196 }
0197
0198 if (debug)
0199 edm::LogVerbatim("IsoTrack") << "eECALmatrix: with " << idEBEE.size() << " EB+EE hits and "
0200 << "spike flag " << flag;
0201 double etot = (!idEBEE.empty()) ? spr::energyECAL(idEBEE, hitsEB, hitsEE, ebThr, eeThr, tMin, tMax, debug) : 0;
0202 return std::pair<double, bool>(etot, flag);
0203 }
0204
0205 std::pair<double, bool> eECALmatrix(const DetId& detId,
0206 edm::Handle<EcalRecHitCollection>& hitsEB,
0207 edm::Handle<EcalRecHitCollection>& hitsEE,
0208 const EcalChannelStatus& chStatus,
0209 const CaloGeometry* geo,
0210 const CaloTopology* caloTopology,
0211 const EcalSeverityLevelAlgo* sevlv,
0212 const EcalPFRecHitThresholds* eThresholds,
0213 int ieta,
0214 int iphi,
0215 bool debug) {
0216 std::vector<DetId> vdets;
0217 spr::matrixECALIds(detId, ieta, iphi, geo, caloTopology, vdets, debug);
0218 if (debug)
0219 edm::LogVerbatim("IsoTrack") << "Inside eECALmatrix " << 2 * ieta + 1 << "X" << 2 * iphi + 1 << " nXtals "
0220 << vdets.size();
0221 bool flag(true);
0222 for (const auto& id : vdets) {
0223 if ((id.det() == DetId::Ecal) && (id.subdetId() == EcalBarrel)) {
0224 if (sevlv->severityLevel(id, (*hitsEB)) == EcalSeverityLevel::kWeird)
0225 flag = false;
0226 } else if ((id.det() == DetId::Ecal) && (id.subdetId() == EcalEndcap)) {
0227 if (sevlv->severityLevel(id, (*hitsEE)) == EcalSeverityLevel::kWeird)
0228 flag = false;
0229 }
0230 }
0231 double energySum = 0.0;
0232 for (unsigned int i1 = 0; i1 < vdets.size(); i1++) {
0233 if (vdets[i1] != DetId(0)) {
0234 std::vector<EcalRecHitCollection::const_iterator> hit;
0235 if (vdets[i1].subdetId() == EcalBarrel) {
0236 spr::findHit(hitsEB, vdets[i1], hit, debug);
0237 } else if (vdets[i1].subdetId() == EcalEndcap) {
0238 spr::findHit(hitsEE, vdets[i1], hit, debug);
0239 }
0240 std::ostringstream st1;
0241 if (debug)
0242 st1 << "Crystal 0x" << std::hex << vdets[i1]() << std::dec;
0243 double ener = 0, ethr = static_cast<double>((*eThresholds)[vdets[i1]]);
0244 for (unsigned int ihit = 0; ihit < hit.size(); ihit++) {
0245 double en = 0;
0246 if (vdets[i1].subdetId() == EcalBarrel) {
0247 if (hit[ihit] != hitsEB->end())
0248 en = hit[ihit]->energy();
0249 } else if (vdets[i1].subdetId() == EcalEndcap) {
0250 if (hit[ihit] != hitsEE->end())
0251 en = hit[ihit]->energy();
0252 }
0253 if (debug)
0254 st1 << " " << ihit << " " << en << " Thr " << ethr;
0255 ener += en;
0256 }
0257 if (debug)
0258 edm::LogVerbatim("IsoTrack") << st1.str();
0259 if (ener > ethr)
0260 energySum += ener;
0261 }
0262 }
0263 return std::pair<double, bool>(energySum, flag);
0264 }
0265
0266 }
