Back to home page

Project CMSSW displayed by LXR

 
 

    


File indexing completed on 2021-08-12 03:21:14

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 }  // namespace spr