File indexing completed on 2024-04-06 11:59:21
0001 #include "DataFormats/HcalDetId/interface/HcalDetId.h"
0002
0003 #include "Calibration/IsolatedParticles/interface/MatrixHCALDetIds.h"
0004 #include "Calibration/IsolatedParticles/interface/FindDistCone.h"
0005 #include "Calibration/IsolatedParticles/interface/DebugInfo.h"
0006 #include "Geometry/HcalTowerAlgo/interface/HcalGeometry.h"
0007 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0008
0009 #include <algorithm>
0010
0011 namespace spr {
0012
0013 std::vector<DetId> matrixHCALIds(
0014 std::vector<DetId>& dets, const HcalTopology* topology, int ieta, int iphi, bool includeHO, bool debug) {
0015 if (debug) {
0016 edm::LogVerbatim("IsoTrack") << "matrixHCALIds::Add " << ieta << " rows and " << iphi << " columns of cells for "
0017 << dets.size() << " cells";
0018 spr::debugHcalDets(0, dets);
0019 }
0020
0021 std::vector<DetId> vdetN = spr::newHCALIdNS(dets, 0, topology, true, ieta, iphi, debug);
0022 std::vector<DetId> vdetS = spr::newHCALIdNS(dets, 0, topology, false, ieta, iphi, debug);
0023 for (unsigned int i1 = 0; i1 < vdetS.size(); i1++) {
0024 if (std::count(vdetN.begin(), vdetN.end(), vdetS[i1]) == 0)
0025 vdetN.push_back(vdetS[i1]);
0026 }
0027
0028 vdetS = spr::matrixHCALIdsDepth(vdetN, topology, includeHO, debug);
0029
0030 if (debug) {
0031 edm::LogVerbatim("IsoTrack") << "matrixHCALIds::Total number of cells found is " << vdetS.size();
0032 spr::debugHcalDets(0, vdetS);
0033 }
0034 return vdetS;
0035 }
0036
0037 std::vector<DetId> matrixHCALIds(const DetId& det,
0038 const CaloGeometry* geo,
0039 const HcalTopology* topology,
0040 double dR,
0041 const GlobalVector& trackMom,
0042 bool includeHO,
0043 bool debug) {
0044 HcalDetId hcdet = HcalDetId(det);
0045 GlobalPoint core = (static_cast<const HcalGeometry*>(geo->getSubdetectorGeometry(hcdet)))->getPosition(hcdet);
0046 std::vector<DetId> dets, vdetx;
0047 dets.push_back(det);
0048 int ietaphi = (int)(dR / 15.0) + 1;
0049 std::vector<DetId> vdets = spr::matrixHCALIds(dets, topology, ietaphi, ietaphi, includeHO, debug);
0050 for (unsigned int i = 0; i < vdets.size(); ++i) {
0051 HcalDetId hcdet = HcalDetId(vdets[i]);
0052 GlobalPoint rpoint = (static_cast<const HcalGeometry*>(geo->getSubdetectorGeometry(hcdet)))->getPosition(hcdet);
0053 if (spr::getDistInPlaneTrackDir(core, trackMom, rpoint) < dR) {
0054 vdetx.push_back(vdets[i]);
0055 }
0056 }
0057
0058 if (debug) {
0059 edm::LogVerbatim("IsoTrack") << "matrixHCALIds::Final List of cells for dR " << dR << " is with " << vdetx.size()
0060 << " from original list of " << vdets.size() << " cells";
0061 spr::debugHcalDets(0, vdetx);
0062 }
0063 return vdetx;
0064 }
0065
0066 std::vector<DetId> matrixHCALIds(std::vector<DetId>& dets,
0067 const HcalTopology* topology,
0068 int ietaE,
0069 int ietaW,
0070 int iphiN,
0071 int iphiS,
0072 bool includeHO,
0073 bool debug) {
0074 if (debug) {
0075 edm::LogVerbatim("IsoTrack") << "matrixHCALIds::Add " << ietaE << "|" << ietaW << " rows and " << iphiN << "|"
0076 << iphiS << " columns of cells for " << dets.size() << " cells";
0077 spr::debugHcalDets(0, dets);
0078 }
0079
0080 std::vector<DetId> vdetN = spr::newHCALIdNS(dets, 0, topology, true, ietaE, ietaW, iphiN, iphiS, debug);
0081 std::vector<DetId> vdetS = spr::newHCALIdNS(dets, 0, topology, false, ietaE, ietaW, iphiN, iphiS, debug);
0082 for (unsigned int i1 = 0; i1 < vdetS.size(); i1++) {
0083 if (std::count(vdetN.begin(), vdetN.end(), vdetS[i1]) == 0)
0084 vdetN.push_back(vdetS[i1]);
0085 }
0086
0087 vdetS = spr::matrixHCALIdsDepth(vdetN, topology, includeHO, debug);
0088
0089 if (debug) {
0090 edm::LogVerbatim("IsoTrack") << "matrixHCALIds::Total number of cells found is " << vdetS.size();
0091 spr::debugHcalDets(0, vdetS);
0092 }
0093 return vdetS;
0094 }
0095
0096 std::vector<DetId> newHCALIdNS(std::vector<DetId>& dets,
0097 unsigned int last,
0098 const HcalTopology* topology,
0099 bool shiftNorth,
0100 int ieta,
0101 int iphi,
0102 bool debug) {
0103 if (debug) {
0104 edm::LogVerbatim("IsoTrack") << "newHCALIdNS::Add " << iphi << " columns of cells along " << shiftNorth << " for "
0105 << (dets.size() - last) << " cells";
0106 spr::debugHcalDets(last, dets);
0107 }
0108
0109 std::vector<DetId> vdets;
0110 vdets.insert(vdets.end(), dets.begin(), dets.end());
0111 std::vector<DetId> vdetE, vdetW;
0112 if (last == 0) {
0113 vdetE = spr::newHCALIdEW(dets, last, topology, true, ieta, debug);
0114 vdetW = spr::newHCALIdEW(dets, last, topology, false, ieta, debug);
0115 for (unsigned int i1 = 0; i1 < vdetW.size(); i1++) {
0116 if (std::count(vdets.begin(), vdets.end(), vdetW[i1]) == 0)
0117 vdets.push_back(vdetW[i1]);
0118 }
0119 for (unsigned int i1 = 0; i1 < vdetE.size(); i1++) {
0120 if (std::count(vdets.begin(), vdets.end(), vdetE[i1]) == 0)
0121 vdets.push_back(vdetE[i1]);
0122 }
0123
0124 if (debug) {
0125 edm::LogVerbatim("IsoTrack") << "newHCALIdNS::With Added cells along E/W results a set of "
0126 << (vdets.size() - dets.size()) << " new cells";
0127 spr::debugHcalDets(dets.size(), vdets);
0128 }
0129 }
0130 unsigned int last0 = vdets.size();
0131 if (iphi > 0) {
0132 std::vector<DetId> vdetnew;
0133 for (unsigned int i1 = last; i1 < dets.size(); i1++) {
0134 std::vector<DetId> vdet;
0135 if (shiftNorth)
0136 vdet = topology->north(dets[i1]);
0137 else
0138 vdet = topology->south(dets[i1]);
0139 for (unsigned int i2 = 0; i2 < vdet.size(); i2++) {
0140 if (std::count(vdets.begin(), vdets.end(), vdet[i2]) == 0)
0141 vdetnew.push_back(vdet[i2]);
0142 }
0143 }
0144 iphi--;
0145 vdetE = spr::newHCALIdEW(vdetnew, 0, topology, true, ieta, debug);
0146 vdetW = spr::newHCALIdEW(vdetnew, 0, topology, false, ieta, debug);
0147 for (unsigned int i2 = 0; i2 < vdetW.size(); i2++) {
0148 if (std::count(vdets.begin(), vdets.end(), vdetW[i2]) == 0 &&
0149 std::count(vdetnew.begin(), vdetnew.end(), vdetW[i2]) == 0)
0150 vdets.push_back(vdetW[i2]);
0151 }
0152 for (unsigned int i2 = 0; i2 < vdetE.size(); i2++) {
0153 if (std::count(vdets.begin(), vdets.end(), vdetE[i2]) == 0 &&
0154 std::count(vdetnew.begin(), vdetnew.end(), vdetE[i2]) == 0)
0155 vdets.push_back(vdetE[i2]);
0156 }
0157 last = vdets.size();
0158 vdets.insert(vdets.end(), vdetnew.begin(), vdetnew.end());
0159
0160 if (debug) {
0161 edm::LogVerbatim("IsoTrack") << "newHCALIdNS::Addition results a set of " << (vdets.size() - last0)
0162 << " new cells";
0163 spr::debugHcalDets(last0, vdets);
0164 }
0165 last0 = last;
0166 }
0167
0168 if (iphi > 0) {
0169 last = last0;
0170 return spr::newHCALIdNS(vdets, last, topology, shiftNorth, ieta, iphi, debug);
0171 } else {
0172 if (debug) {
0173 edm::LogVerbatim("IsoTrack") << "newHCALIdNS::Final list consists of " << vdets.size() << " cells";
0174 spr::debugHcalDets(0, vdets);
0175 }
0176 return vdets;
0177 }
0178 }
0179
0180 std::vector<DetId> newHCALIdNS(std::vector<DetId>& dets,
0181 unsigned int last,
0182 const HcalTopology* topology,
0183 bool shiftNorth,
0184 int ietaE,
0185 int ietaW,
0186 int iphiN,
0187 int iphiS,
0188 bool debug) {
0189 if (debug) {
0190 edm::LogVerbatim("IsoTrack") << "newHCALIdNS::Add " << iphiN << "|" << iphiS << " columns of cells along "
0191 << shiftNorth << " for " << (dets.size() - last) << " cells";
0192 spr::debugHcalDets(last, dets);
0193 }
0194
0195 std::vector<DetId> vdets;
0196 vdets.insert(vdets.end(), dets.begin(), dets.end());
0197 std::vector<DetId> vdetE, vdetW;
0198 if (last == 0) {
0199 vdetE = spr::newHCALIdEW(dets, last, topology, true, ietaE, ietaW, debug);
0200 vdetW = spr::newHCALIdEW(dets, last, topology, false, ietaE, ietaW, debug);
0201 for (unsigned int i1 = 0; i1 < vdetW.size(); i1++) {
0202 if (std::count(vdets.begin(), vdets.end(), vdetW[i1]) == 0)
0203 vdets.push_back(vdetW[i1]);
0204 }
0205 for (unsigned int i1 = 0; i1 < vdetE.size(); i1++) {
0206 if (std::count(vdets.begin(), vdets.end(), vdetE[i1]) == 0)
0207 vdets.push_back(vdetE[i1]);
0208 }
0209
0210 if (debug) {
0211 edm::LogVerbatim("IsoTrack") << "newHCALIdNS::With Added cells along E/W results a set of "
0212 << (vdets.size() - dets.size()) << " new cells";
0213 spr::debugHcalDets(dets.size(), vdets);
0214 }
0215 }
0216 unsigned int last0 = vdets.size();
0217 int iphi = iphiS;
0218 if (shiftNorth)
0219 iphi = iphiN;
0220 if (iphi > 0) {
0221 std::vector<DetId> vdetnew;
0222 for (unsigned int i1 = last; i1 < dets.size(); i1++) {
0223 std::vector<DetId> vdet;
0224 if (shiftNorth)
0225 vdet = topology->north(dets[i1]);
0226 else
0227 vdet = topology->south(dets[i1]);
0228 for (unsigned int i2 = 0; i2 < vdet.size(); i2++) {
0229 if (std::count(vdets.begin(), vdets.end(), vdet[i2]) == 0)
0230 vdetnew.push_back(vdet[i2]);
0231 }
0232 }
0233 iphi--;
0234 vdetE = spr::newHCALIdEW(vdetnew, 0, topology, true, ietaE, ietaW, debug);
0235 vdetW = spr::newHCALIdEW(vdetnew, 0, topology, false, ietaE, ietaW, debug);
0236 for (unsigned int i2 = 0; i2 < vdetW.size(); i2++) {
0237 if (std::count(vdets.begin(), vdets.end(), vdetW[i2]) == 0 &&
0238 std::count(vdetnew.begin(), vdetnew.end(), vdetW[i2]) == 0)
0239 vdets.push_back(vdetW[i2]);
0240 }
0241 for (unsigned int i2 = 0; i2 < vdetE.size(); i2++) {
0242 if (std::count(vdets.begin(), vdets.end(), vdetE[i2]) == 0 &&
0243 std::count(vdetnew.begin(), vdetnew.end(), vdetE[i2]) == 0)
0244 vdets.push_back(vdetE[i2]);
0245 }
0246 last = vdets.size();
0247 vdets.insert(vdets.end(), vdetnew.begin(), vdetnew.end());
0248
0249 if (debug) {
0250 edm::LogVerbatim("IsoTrack") << "newHCALIdNS::Addition results a set of " << (vdets.size() - last0)
0251 << " new cells";
0252 spr::debugHcalDets(last0, vdets);
0253 }
0254 last0 = last;
0255 }
0256 if (shiftNorth)
0257 iphiN = iphi;
0258 else
0259 iphiS = iphi;
0260
0261 if (iphi > 0) {
0262 last = last0;
0263 return spr::newHCALIdNS(vdets, last, topology, shiftNorth, ietaE, ietaW, iphiN, iphiS, debug);
0264 } else {
0265 if (debug) {
0266 edm::LogVerbatim("IsoTrack") << "newHCALIdNS::Final list consists of " << vdets.size() << " cells";
0267 spr::debugHcalDets(0, vdets);
0268 }
0269 return vdets;
0270 }
0271 }
0272
0273 std::vector<DetId> newHCALIdEW(
0274 std::vector<DetId>& dets, unsigned int last, const HcalTopology* topology, bool shiftEast, int ieta, bool debug) {
0275 if (debug) {
0276 edm::LogVerbatim("IsoTrack") << "newHCALIdEW::Add " << ieta << " rows of cells along " << shiftEast << " for "
0277 << (dets.size() - last) << " cells";
0278 spr::debugHcalDets(last, dets);
0279 }
0280
0281 std::vector<DetId> vdets;
0282 vdets.insert(vdets.end(), dets.begin(), dets.end());
0283 if (ieta > 0) {
0284 for (unsigned int i1 = last; i1 < dets.size(); i1++) {
0285 std::vector<DetId> vdet;
0286 if (shiftEast)
0287 vdet = topology->east(dets[i1]);
0288 else
0289 vdet = topology->west(dets[i1]);
0290 for (unsigned int i2 = 0; i2 < vdet.size(); i2++) {
0291 if (std::count(vdets.begin(), vdets.end(), vdet[i2]) == 0)
0292 vdets.push_back(vdet[i2]);
0293 }
0294 }
0295 ieta--;
0296 }
0297
0298 if (debug) {
0299 edm::LogVerbatim("IsoTrack") << "newHCALIdEW::Addition results a set of " << (vdets.size() - dets.size())
0300 << " new cells";
0301 spr::debugHcalDets(dets.size(), vdets);
0302 }
0303 if (ieta > 0) {
0304 last = dets.size();
0305 return spr::newHCALIdEW(vdets, last, topology, shiftEast, ieta, debug);
0306 } else {
0307 if (debug) {
0308 edm::LogVerbatim("IsoTrack") << "newHCALIdEW::Final list (EW) consists of " << vdets.size() << " cells";
0309 spr::debugHcalDets(0, vdets);
0310 }
0311 return vdets;
0312 }
0313 }
0314
0315 std::vector<DetId> newHCALIdEW(std::vector<DetId>& dets,
0316 unsigned int last,
0317 const HcalTopology* topology,
0318 bool shiftEast,
0319 int ietaE,
0320 int ietaW,
0321 bool debug) {
0322 if (debug) {
0323 edm::LogVerbatim("IsoTrack") << "newHCALIdEW::Add " << ietaE << "|" << ietaW << " rows of cells along "
0324 << shiftEast << " for " << (dets.size() - last) << " cells";
0325 spr::debugHcalDets(last, dets);
0326 }
0327
0328 int ieta = ietaW;
0329 if (shiftEast)
0330 ieta = ietaE;
0331 std::vector<DetId> vdets;
0332 vdets.insert(vdets.end(), dets.begin(), dets.end());
0333 if (ieta > 0) {
0334 for (unsigned int i1 = last; i1 < dets.size(); i1++) {
0335 std::vector<DetId> vdet;
0336 if (shiftEast)
0337 vdet = topology->east(dets[i1]);
0338 else
0339 vdet = topology->west(dets[i1]);
0340 for (unsigned int i2 = 0; i2 < vdet.size(); i2++) {
0341 if (std::count(vdets.begin(), vdets.end(), vdet[i2]) == 0)
0342 vdets.push_back(vdet[i2]);
0343 }
0344 }
0345 ieta--;
0346 }
0347 if (shiftEast)
0348 ietaE = ieta;
0349 else
0350 ietaW = ieta;
0351
0352 if (debug) {
0353 edm::LogVerbatim("IsoTrack") << "newHCALIdEW::Addition results a set of " << (vdets.size() - dets.size())
0354 << " new cells";
0355 spr::debugHcalDets(dets.size(), vdets);
0356 }
0357
0358 if (ieta > 0) {
0359 last = dets.size();
0360 return spr::newHCALIdEW(vdets, last, topology, shiftEast, ietaE, ietaW, debug);
0361 } else {
0362 if (debug) {
0363 edm::LogVerbatim("IsoTrack") << "newHCALIdEW::Final list (EW) consists of " << vdets.size() << " cells";
0364 spr::debugHcalDets(0, vdets);
0365 }
0366 return vdets;
0367 }
0368 }
0369
0370 std::vector<DetId> matrixHCALIdsDepth(std::vector<DetId>& dets,
0371 const HcalTopology* topology,
0372 bool includeHO,
0373 bool debug) {
0374 if (debug) {
0375 edm::LogVerbatim("IsoTrack") << "matrixHCALIdsDepth::Add cells with higher depths with HO"
0376 << "Flag set to " << includeHO << " to existing " << dets.size() << " cells";
0377 spr::debugHcalDets(0, dets);
0378 }
0379 std::vector<DetId> vdets(dets);
0380 for (unsigned int i1 = 0; i1 < dets.size(); i1++) {
0381 HcalDetId vdet = dets[i1];
0382 for (int idepth = 0; idepth < 3; idepth++) {
0383 std::vector<DetId> vUpDetId = topology->up(vdet);
0384 if (!vUpDetId.empty()) {
0385 if (includeHO || vUpDetId[0].subdetId() != (int)(HcalOuter)) {
0386 int n = std::count(vdets.begin(), vdets.end(), vUpDetId[0]);
0387 if (n == 0) {
0388 if (debug)
0389 edm::LogVerbatim("IsoTrack")
0390 << "matrixHCALIdsDepth:: Depth " << idepth << " " << vdet << " " << (HcalDetId)vUpDetId[0];
0391 vdets.push_back(vUpDetId[0]);
0392 }
0393 }
0394 vdet = vUpDetId[0];
0395 }
0396 }
0397 }
0398
0399 if (debug) {
0400 edm::LogVerbatim("IsoTrack") << "matrixHCALIdsDepth::Final list contains " << vdets.size() << " cells";
0401 spr::debugHcalDets(0, vdets);
0402 }
0403 return vdets;
0404 }
0405
0406 }
