File indexing completed on 2024-04-06 11:59:20
0001 #include "Calibration/IsolatedParticles/interface/DebugInfo.h"
0002 #include "Calibration/IsolatedParticles/interface/eHCALMatrix.h"
0003 #include "DataFormats/HcalDetId/interface/HcalSubdetector.h"
0004 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0005
0006 #include <algorithm>
0007 #include <sstream>
0008
0009 namespace spr {
0010 double eHCALmatrix(const HcalTopology* topology,
0011 const DetId& det0,
0012 std::vector<PCaloHit>& hits,
0013 int ieta,
0014 int iphi,
0015 bool includeHO,
0016 double hbThr,
0017 double heThr,
0018 double hfThr,
0019 double hoThr,
0020 double tMin,
0021 double tMax,
0022 bool debug) {
0023 HcalDetId hcid0(det0.rawId());
0024 HcalDetId hcid(hcid0.subdet(), hcid0.ieta(), hcid0.iphi(), 1);
0025 DetId det(hcid.rawId());
0026 if (debug)
0027 edm::LogVerbatim("IsoTrack") << "Inside eHCALmatrix " << 2 * ieta + 1 << "X" << 2 * iphi + 1
0028 << " Inclusion of HO Flag " << includeHO;
0029
0030 double energySum(0);
0031 std::vector<DetId> dets(1, det);
0032 std::vector<DetId> vdets = spr::matrixHCALIds(dets, topology, ieta, iphi, includeHO, false);
0033 if (debug) {
0034 edm::LogVerbatim("IsoTrack") << "matrixHCALIds::Total number of cells found is " << vdets.size();
0035 spr::debugHcalDets(0, vdets);
0036 }
0037
0038 int khit(0);
0039 for (unsigned int i = 0; i < vdets.size(); i++) {
0040 std::vector<std::vector<PCaloHit>::const_iterator> hit = spr::findHit(hits, vdets[i]);
0041 double energy = 0;
0042 int subdet = ((HcalDetId)(vdets[i].rawId())).subdet();
0043 double eThr = spr::eHCALThreshold(subdet, hbThr, heThr, hfThr, hoThr);
0044 for (unsigned int ihit = 0; ihit < hit.size(); ihit++) {
0045 if (hit[ihit] != hits.end()) {
0046 khit++;
0047 if (debug)
0048 edm::LogVerbatim("IsoTrack") << "energyHCAL:: Hit " << khit << " " << (HcalDetId)vdets[i] << " E "
0049 << hit[ihit]->energy() << " t " << hit[ihit]->time();
0050
0051 if (hit[ihit]->time() > tMin && hit[ihit]->time() < tMax) {
0052 energy += hit[ihit]->energy();
0053 }
0054 }
0055 }
0056 if (energy > eThr)
0057 energySum += energy;
0058 }
0059
0060 if (debug)
0061 edm::LogVerbatim("IsoTrack") << "eHCALmatrix::Total energy " << energySum;
0062 return energySum;
0063 }
0064
0065 double eHCALmatrix(const CaloGeometry* geo,
0066 const HcalTopology* topology,
0067 const DetId& det0,
0068 std::vector<PCaloHit>& hits,
0069 int ieta,
0070 int iphi,
0071 HcalDetId& hotCell,
0072 bool includeHO,
0073 bool debug) {
0074 HcalDetId hcid0(det0.rawId());
0075 HcalDetId hcid(hcid0.subdet(), hcid0.ieta(), hcid0.iphi(), 1);
0076 DetId det(hcid.rawId());
0077 std::vector<DetId> dets(1, det);
0078 std::vector<DetId> vdets = spr::matrixHCALIds(dets, topology, ieta, iphi, includeHO, debug);
0079 hotCell = hcid0;
0080
0081 std::vector<std::vector<PCaloHit>::const_iterator> hitlist;
0082 for (unsigned int i = 0; i < vdets.size(); i++) {
0083 std::vector<std::vector<PCaloHit>::const_iterator> hit = spr::findHit(hits, vdets[i]);
0084 hitlist.insert(hitlist.end(), hit.begin(), hit.end());
0085 }
0086
0087 double energySum(0);
0088 for (unsigned int ihit = 0; ihit < hitlist.size(); ihit++)
0089 energySum += hitlist[ihit]->energy();
0090
0091
0092 dets.clear();
0093 std::vector<double> energies;
0094 for (unsigned int ihit = 0; ihit < hitlist.size(); ihit++) {
0095 double energy = hitlist[ihit]->energy();
0096 HcalDetId id0 = HcalDetId(hitlist[ihit]->id());
0097 if ((id0.subdet() != HcalOuter) || includeHO) {
0098 HcalDetId id1(id0.subdet(), id0.ieta(), id0.iphi(), 1);
0099 bool found(false);
0100 for (unsigned int idet = 0; idet < dets.size(); ++idet) {
0101 if (id1 == HcalDetId(dets[idet])) {
0102 energies[idet] += energy;
0103 found = true;
0104 break;
0105 }
0106 }
0107 if (!found) {
0108 dets.push_back(DetId(id1));
0109 energies.push_back(energy);
0110 }
0111 }
0112 }
0113 double energyMax(-99.);
0114 for (unsigned int ihit = 0; ihit < dets.size(); ihit++) {
0115 if (energies[ihit] > energyMax) {
0116 energyMax = energies[ihit];
0117 hotCell = HcalDetId(dets[ihit]);
0118 }
0119 }
0120 return energySum;
0121 }
0122
0123 void energyHCALCell(HcalDetId detId,
0124 std::vector<PCaloHit>& hits,
0125 std::vector<std::pair<double, int> >& energyCell,
0126 int maxDepth,
0127 double hbThr,
0128 double heThr,
0129 double hfThr,
0130 double hoThr,
0131 double tMin,
0132 double tMax,
0133 int depthHE,
0134 bool debug) {
0135 energyCell.clear();
0136 int subdet = detId.subdet();
0137 double eThr = spr::eHCALThreshold(subdet, hbThr, heThr, hfThr, hoThr);
0138 bool hbhe = (detId.ietaAbs() == 16);
0139 if (debug)
0140 edm::LogVerbatim("IsoTrack") << "energyHCALCell: input ID " << detId << " MaxDepth " << maxDepth
0141 << " Threshold (E) " << eThr << " (T) " << tMin << ":" << tMax;
0142
0143 for (int i = 0; i < maxDepth; i++) {
0144 HcalSubdetector subdet0 = (hbhe) ? ((i + 1 >= depthHE) ? HcalEndcap : HcalBarrel) : detId.subdet();
0145 HcalDetId hcid(subdet0, detId.ieta(), detId.iphi(), i + 1);
0146 DetId det(hcid.rawId());
0147 std::vector<std::vector<PCaloHit>::const_iterator> hit = spr::findHit(hits, det);
0148 double energy(0);
0149 for (unsigned int ihit = 0; ihit < hit.size(); ++ihit) {
0150 if (hit[ihit]->time() > tMin && hit[ihit]->time() < tMax)
0151 energy += hit[ihit]->energy();
0152 if (debug)
0153 edm::LogVerbatim("IsoTrack") << "energyHCALCell:: Hit[" << ihit << "] " << hcid << " E "
0154 << hit[ihit]->energy() << " t " << hit[ihit]->time();
0155 }
0156
0157 if (debug)
0158 edm::LogVerbatim("IsoTrack") << "energyHCALCell:: Cell " << hcid << " E " << energy << " from " << hit.size()
0159 << " threshold " << eThr;
0160 if (energy > eThr && !hit.empty()) {
0161 energyCell.push_back(std::pair<double, int>(energy, i + 1));
0162 }
0163 }
0164
0165 if (debug) {
0166 std::ostringstream st1;
0167 st1 << "energyHCALCell:: " << energyCell.size() << " entries from " << maxDepth << " depths:";
0168 for (unsigned int i = 0; i < energyCell.size(); ++i) {
0169 st1 << " [" << i << "] (" << energyCell[i].first << ":" << energyCell[i].second << ")";
0170 }
0171 edm::LogVerbatim("IsoTrack") << st1.str();
0172 }
0173 }
0174
0175 HcalDetId getHotCell(std::vector<HBHERecHitCollection::const_iterator>& hit, bool includeHO, int useRaw, bool) {
0176 std::vector<HcalDetId> dets;
0177 std::vector<double> energies;
0178 for (unsigned int ihit = 0; ihit < hit.size(); ihit++) {
0179 double energy = getRawEnergy(hit.at(ihit), useRaw);
0180 HcalDetId id0 = hit.at(ihit)->id();
0181 if ((id0.subdet() != HcalOuter) || includeHO) {
0182 HcalDetId id1(id0.subdet(), id0.ieta(), id0.iphi(), 1);
0183 bool found(false);
0184 for (unsigned int idet = 0; idet < dets.size(); ++idet) {
0185 if (id1 == dets[idet]) {
0186 energies[idet] += energy;
0187 found = true;
0188 break;
0189 }
0190 }
0191 if (!found) {
0192 dets.push_back(id1);
0193 energies.push_back(energy);
0194 }
0195 }
0196 }
0197 double energyMax(-99.);
0198 HcalDetId hotCell;
0199 for (unsigned int ihit = 0; ihit < dets.size(); ihit++) {
0200 if (energies[ihit] > energyMax) {
0201 energyMax = energies[ihit];
0202 hotCell = dets[ihit];
0203 }
0204 }
0205 return hotCell;
0206 }
0207
0208 HcalDetId getHotCell(std::vector<std::vector<PCaloHit>::const_iterator>& hit, bool includeHO, int useRaw, bool) {
0209 std::vector<HcalDetId> dets;
0210 std::vector<double> energies;
0211 for (unsigned int ihit = 0; ihit < hit.size(); ihit++) {
0212 double energy = hit.at(ihit)->energy();
0213 HcalDetId id0 = getRawEnergy(hit.at(ihit), useRaw);
0214 if ((id0.subdet() != HcalOuter) || includeHO) {
0215 HcalDetId id1(id0.subdet(), id0.ieta(), id0.iphi(), 1);
0216 bool found(false);
0217 for (unsigned int idet = 0; idet < dets.size(); ++idet) {
0218 if (id1 == dets[idet]) {
0219 energies[idet] += energy;
0220 found = true;
0221 break;
0222 }
0223 }
0224 if (!found) {
0225 dets.push_back(id1);
0226 energies.push_back(energy);
0227 }
0228 }
0229 }
0230 double energyMax(-99.);
0231 HcalDetId hotCell;
0232 for (unsigned int ihit = 0; ihit < dets.size(); ihit++) {
0233 if (energies[ihit] > energyMax) {
0234 energyMax = energies[ihit];
0235 hotCell = dets[ihit];
0236 }
0237 }
0238 return hotCell;
0239 }
0240
0241 double eHCALThreshold(int subdet, double hbThr, double heThr, double hfThr, double hoThr) {
0242 double eThr = hbThr;
0243 if (subdet == (int)(HcalEndcap))
0244 eThr = heThr;
0245 else if (subdet == (int)(HcalForward))
0246 eThr = hfThr;
0247 else if (subdet == (int)(HcalOuter))
0248 eThr = hoThr;
0249 return eThr;
0250 }
0251 }
