Back to home page

Project CMSSW displayed by LXR
 
 

    

File indexing completed on 2024-04-06 11:59:20 

0001 #include "DataFormats/EcalDetId/interface/EBDetId.h"
0002 #include "DataFormats/EcalDetId/interface/EcalSubdetector.h"
0003 #include "DataFormats/EcalDetId/interface/EEDetId.h"
0004 
0005 #include "Calibration/IsolatedParticles/interface/MatrixECALDetIds.h"
0006 #include "Calibration/IsolatedParticles/interface/FindDistCone.h"
0007 #include "Calibration/IsolatedParticles/interface/DebugInfo.h"
0008 
0009 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0010 
0011 #include <algorithm>
0012 #include <iostream>
0013 
0014 namespace spr {
0015 
0016   void matrixECALIds(const DetId& det,
0017                      int ieta,
0018                      int iphi,
0019                      const CaloGeometry* geo,
0020                      const CaloTopology* caloTopology,
0021                      std::vector<DetId>& vdets,
0022                      bool debug,
0023                      bool ignoreTransition) {
0024     const CaloSubdetectorTopology* barrelTopo = (caloTopology->getSubdetectorTopology(DetId::Ecal, EcalBarrel));
0025     const CaloSubdetectorTopology* endcapTopo = (caloTopology->getSubdetectorTopology(DetId::Ecal, EcalEndcap));
0026     const EcalBarrelGeometry* barrelGeom =
0027         (dynamic_cast<const EcalBarrelGeometry*>(geo->getSubdetectorGeometry(DetId::Ecal, EcalBarrel)));
0028     const EcalEndcapGeometry* endcapGeom =
0029         (dynamic_cast<const EcalEndcapGeometry*>(geo->getSubdetectorGeometry(DetId::Ecal, EcalEndcap)));
0030 
0031     if (debug) {
0032       edm::LogVerbatim("IsoTrack") << "matrixECALIds::Add " << ieta << " rows and " << iphi
0033                                    << " columns of cells for 1 cell" << spr::debugEcalDet(0, det).str();
0034     }
0035     std::vector<DetId> dets(1, det);
0036     std::vector<CaloDirection> dirs(1, NORTH);
0037     vdets = spr::newECALIdNS(
0038         dets, 0, ieta, iphi, dirs, barrelTopo, endcapTopo, barrelGeom, endcapGeom, debug, ignoreTransition);
0039     dirs[0] = SOUTH;
0040     std::vector<DetId> vdetS = spr::newECALIdNS(
0041         dets, 0, ieta, iphi, dirs, barrelTopo, endcapTopo, barrelGeom, endcapGeom, debug, ignoreTransition);
0042     for (unsigned int i1 = 0; i1 < vdetS.size(); i1++) {
0043       if (std::count(vdets.begin(), vdets.end(), vdetS[i1]) == 0)
0044         vdets.push_back(vdetS[i1]);
0045     }
0046     unsigned int ndet = (2 * ieta + 1) * (2 * iphi + 1);
0047     if (vdets.size() != ndet) {
0048       std::vector<DetId> vdetExtra;
0049       spr::extraIds(det, vdets, ieta, ieta, iphi, iphi, barrelGeom, endcapGeom, vdetExtra, debug);
0050       if (!vdetExtra.empty())
0051         vdets.insert(vdets.end(), vdetExtra.begin(), vdetExtra.end());
0052     }
0053 
0054     if (debug) {
0055       edm::LogVerbatim("IsoTrack") << "matrixECALIds::Total number of cells found is " << vdets.size();
0056       spr::debugEcalDets(0, vdets);
0057     }
0058   }
0059 
0060   std::vector<DetId> matrixECALIds(const DetId& det,
0061                                    int ieta,
0062                                    int iphi,
0063                                    const CaloGeometry* geo,
0064                                    const CaloTopology* caloTopology,
0065                                    bool debug,
0066                                    bool ignoreTransition) {
0067     std::vector<DetId> vdets;
0068     spr::matrixECALIds(det, ieta, iphi, geo, caloTopology, vdets, debug, ignoreTransition);
0069     return vdets;
0070   }
0071 
0072   std::vector<DetId> matrixECALIds(const DetId& det,
0073                                    double dR,
0074                                    const GlobalVector& trackMom,
0075                                    const CaloGeometry* geo,
0076                                    const CaloTopology* caloTopology,
0077                                    bool debug,
0078                                    bool ignoreTransition) {
0079     GlobalPoint core;
0080     if (det.subdetId() == EcalEndcap) {
0081       EEDetId EEid = EEDetId(det);
0082       core = geo->getPosition(EEid);
0083     } else {
0084       EBDetId EBid = EBDetId(det);
0085       core = geo->getPosition(EBid);
0086     }
0087     int ietaphi = (int)(dR / 2.0) + 1;
0088     std::vector<DetId> vdets, vdetx;
0089     spr::matrixECALIds(det, ietaphi, ietaphi, geo, caloTopology, vdets, debug, ignoreTransition);
0090     for (unsigned int i = 0; i < vdets.size(); ++i) {
0091       GlobalPoint rpoint;
0092       if (vdets[i].subdetId() == EcalEndcap) {
0093         EEDetId EEid = EEDetId(vdets[i]);
0094         rpoint = geo->getPosition(EEid);
0095       } else {
0096         EBDetId EBid = EBDetId(vdets[i]);
0097         rpoint = geo->getPosition(EBid);
0098       }
0099       if (spr::getDistInPlaneTrackDir(core, trackMom, rpoint) < dR) {
0100         vdetx.push_back(vdets[i]);
0101       }
0102     }
0103 
0104     if (debug) {
0105       edm::LogVerbatim("IsoTrack") << "matrixECALIds::Final List of cells for dR " << dR << " is with " << vdetx.size()
0106                                    << " from original list of " << vdets.size();
0107       spr::debugEcalDets(0, vdetx);
0108     }
0109     return vdetx;
0110   }
0111 
0112   void matrixECALIds(const DetId& det,
0113                      int ietaE,
0114                      int ietaW,
0115                      int iphiN,
0116                      int iphiS,
0117                      const CaloGeometry* geo,
0118                      const CaloTopology* caloTopology,
0119                      std::vector<DetId>& vdets,
0120                      bool debug,
0121                      bool ignoreTransition) {
0122     const CaloSubdetectorTopology* barrelTopo = (caloTopology->getSubdetectorTopology(DetId::Ecal, EcalBarrel));
0123     const CaloSubdetectorTopology* endcapTopo = (caloTopology->getSubdetectorTopology(DetId::Ecal, EcalEndcap));
0124     const EcalBarrelGeometry* barrelGeom =
0125         (dynamic_cast<const EcalBarrelGeometry*>(geo->getSubdetectorGeometry(DetId::Ecal, EcalBarrel)));
0126     const EcalEndcapGeometry* endcapGeom =
0127         (dynamic_cast<const EcalEndcapGeometry*>(geo->getSubdetectorGeometry(DetId::Ecal, EcalEndcap)));
0128 
0129     if (debug) {
0130       edm::LogVerbatim("IsoTrack") << "matrixECALIds::Add " << ietaE << "|" << ietaW << " rows and " << iphiN << "|"
0131                                    << iphiS << " columns of cells for 1 cell" << spr::debugEcalDet(0, det).str();
0132     }
0133     std::vector<DetId> dets(1, det);
0134     std::vector<CaloDirection> dirs(1, NORTH);
0135     std::vector<int> jetaE(1, ietaE), jetaW(1, ietaW);
0136     std::vector<int> jphiN(1, iphiN), jphiS(1, iphiS);
0137     vdets = spr::newECALIdNS(dets,
0138                              0,
0139                              jetaE,
0140                              jetaW,
0141                              jphiN,
0142                              jphiS,
0143                              dirs,
0144                              barrelTopo,
0145                              endcapTopo,
0146                              barrelGeom,
0147                              endcapGeom,
0148                              debug,
0149                              ignoreTransition);
0150     dirs[0] = SOUTH;
0151     std::vector<DetId> vdetS = spr::newECALIdNS(dets,
0152                                                 0,
0153                                                 jetaE,
0154                                                 jetaW,
0155                                                 jphiN,
0156                                                 jphiS,
0157                                                 dirs,
0158                                                 barrelTopo,
0159                                                 endcapTopo,
0160                                                 barrelGeom,
0161                                                 endcapGeom,
0162                                                 debug,
0163                                                 ignoreTransition);
0164     for (unsigned int i1 = 0; i1 < vdetS.size(); i1++) {
0165       if (std::count(vdets.begin(), vdets.end(), vdetS[i1]) == 0)
0166         vdets.push_back(vdetS[i1]);
0167     }
0168 
0169     unsigned int ndet = (ietaE + ietaW + 1) * (iphiN + iphiS + 1);
0170     if (vdets.size() != ndet) {
0171       std::vector<DetId> vdetExtra;
0172       spr::extraIds(det, vdets, ietaE, ietaW, iphiN, iphiS, barrelGeom, endcapGeom, vdetExtra, debug);
0173       if (!vdetExtra.empty())
0174         vdets.insert(vdets.end(), vdetExtra.begin(), vdetExtra.end());
0175     }
0176 
0177     if (debug) {
0178       edm::LogVerbatim("IsoTrack") << "matrixECALIds::Total number of cells found is " << vdets.size();
0179       spr::debugEcalDets(0, vdets);
0180     }
0181   }
0182 
0183   std::vector<DetId> matrixECALIds(const DetId& det,
0184                                    int ietaE,
0185                                    int ietaW,
0186                                    int iphiN,
0187                                    int iphiS,
0188                                    const CaloGeometry* geo,
0189                                    const CaloTopology* caloTopology,
0190                                    bool debug,
0191                                    bool ignoreTransition) {
0192     std::vector<DetId> vdets;
0193     spr::matrixECALIds(det, ietaE, ietaW, iphiN, iphiS, geo, caloTopology, vdets, debug, ignoreTransition);
0194     return vdets;
0195   }
0196 
0197   std::vector<DetId> newECALIdNS(std::vector<DetId>& dets,
0198                                  unsigned int last,
0199                                  int ieta,
0200                                  int iphi,
0201                                  std::vector<CaloDirection>& dir,
0202                                  const CaloSubdetectorTopology* barrelTopo,
0203                                  const CaloSubdetectorTopology* endcapTopo,
0204                                  const EcalBarrelGeometry* barrelGeom,
0205                                  const EcalEndcapGeometry* endcapGeom,
0206                                  bool debug,
0207                                  bool ignoreTransition) {
0208     if (debug) {
0209       edm::LogVerbatim("IsoTrack") << "newECALIdNS::Add " << iphi << " columns of cells for " << (dets.size() - last)
0210                                    << " cells (last " << last << ")";
0211       spr::debugEcalDets(last, dets, dir);
0212     }
0213     std::vector<DetId> vdets;
0214     std::vector<CaloDirection> dirs;
0215     vdets.insert(vdets.end(), dets.begin(), dets.end());
0216     dirs.insert(dirs.end(), dir.begin(), dir.end());
0217 
0218     std::vector<DetId> vdetE, vdetW;
0219     if (last == 0) {
0220       unsigned int ndet = vdets.size();
0221       std::vector<CaloDirection> dirE(ndet, EAST), dirW(ndet, WEST);
0222       vdetE = spr::newECALIdEW(
0223           dets, last, ieta, dirE, barrelTopo, endcapTopo, barrelGeom, endcapGeom, debug, ignoreTransition);
0224       vdetW = spr::newECALIdEW(
0225           dets, last, ieta, dirW, barrelTopo, endcapTopo, barrelGeom, endcapGeom, debug, ignoreTransition);
0226       for (unsigned int i1 = 0; i1 < vdetW.size(); i1++) {
0227         if (std::count(vdets.begin(), vdets.end(), vdetW[i1]) == 0) {
0228           vdets.push_back(vdetW[i1]);
0229           dirs.push_back(dir[0]);
0230         }
0231       }
0232       for (unsigned int i1 = 0; i1 < vdetE.size(); i1++) {
0233         if (std::count(vdets.begin(), vdets.end(), vdetE[i1]) == 0) {
0234           vdets.push_back(vdetE[i1]);
0235           dirs.push_back(dir[0]);
0236         }
0237       }
0238       if (debug) {
0239         edm::LogVerbatim("IsoTrack") << "newECALIdNS::With Added cells along E/W results a set of "
0240                                      << (vdets.size() - dets.size()) << " new  cells";
0241         spr::debugEcalDets(dets.size(), vdets, dirs);
0242       }
0243     }
0244 
0245     unsigned int last0 = vdets.size();
0246     std::vector<DetId> vdetnew;
0247     std::vector<CaloDirection> dirnew;
0248     if (iphi > 0) {
0249       std::vector<DetId> vdetn(1);
0250       std::vector<CaloDirection> dirn(1);
0251       std::vector<CaloDirection> dirnE(1, EAST), dirnW(1, WEST);
0252       int flag = 0;
0253       for (unsigned int i1 = last; i1 < dets.size(); i1++) {
0254         std::vector<DetId> cells;
0255         spr::simpleMove(
0256             dets[i1], dir[i1], barrelTopo, endcapTopo, barrelGeom, endcapGeom, cells, flag, debug, ignoreTransition);
0257         if (flag != 0) {
0258           if (std::count(vdets.begin(), vdets.end(), cells[0]) == 0) {
0259             vdetn[0] = cells[0];
0260             vdetnew.push_back(vdetn[0]);
0261             dirn[0] = dir[i1];
0262             if (flag < 0) {
0263               if (dirn[0] == NORTH)
0264                 dirn[0] = SOUTH;
0265               else
0266                 dirn[0] = NORTH;
0267             }
0268             dirnew.push_back(dirn[0]);
0269             vdetE = spr::newECALIdEW(
0270                 vdetn, 0, ieta, dirnE, barrelTopo, endcapTopo, barrelGeom, endcapGeom, debug, ignoreTransition);
0271             vdetW = spr::newECALIdEW(
0272                 vdetn, 0, ieta, dirnW, barrelTopo, endcapTopo, barrelGeom, endcapGeom, debug, ignoreTransition);
0273             for (unsigned int i2 = 0; i2 < vdetW.size(); i2++) {
0274               if (std::count(vdets.begin(), vdets.end(), vdetW[i2]) == 0 &&
0275                   std::count(vdetnew.begin(), vdetnew.end(), vdetW[i2]) == 0) {
0276                 vdets.push_back(vdetW[i2]);
0277                 dirs.push_back(dirn[0]);
0278               }
0279             }
0280             for (unsigned int i2 = 0; i2 < vdetE.size(); i2++) {
0281               if (std::count(vdets.begin(), vdets.end(), vdetE[i2]) == 0 &&
0282                   std::count(vdetnew.begin(), vdetnew.end(), vdetE[i2]) == 0) {
0283                 vdets.push_back(vdetE[i2]);
0284                 dirs.push_back(dirn[0]);
0285               }
0286             }
0287           }
0288         }
0289       }
0290       iphi--;
0291       last = vdets.size();
0292       for (unsigned int i2 = 0; i2 < vdetnew.size(); i2++) {
0293         if (std::count(vdets.begin(), vdets.end(), vdetnew[i2]) == 0) {
0294           vdets.push_back(vdetnew[i2]);
0295           dirs.push_back(dirnew[i2]);
0296         }
0297       }
0298       if (debug) {
0299         edm::LogVerbatim("IsoTrack") << "newECALIdNS::Addition results a set of " << (vdets.size() - last0)
0300                                      << " new  cells (last " << last0 << ", iphi " << iphi << ")";
0301         spr::debugEcalDets(last0, vdets, dirs);
0302       }
0303       last0 = last;
0304     }
0305 
0306     if (iphi > 0) {
0307       last = last0;
0308       return spr::newECALIdNS(
0309           vdets, last, ieta, iphi, dirs, barrelTopo, endcapTopo, barrelGeom, endcapGeom, debug, ignoreTransition);
0310     } else {
0311       if (debug) {
0312         edm::LogVerbatim("IsoTrack") << "newECALIdNS::Final list consists of " << vdets.size() << " cells";
0313         spr::debugEcalDets(0, vdets);
0314       }
0315       return vdets;
0316     }
0317   }
0318 
0319   std::vector<DetId> newECALIdNS(std::vector<DetId>& dets,
0320                                  unsigned int last,
0321                                  std::vector<int>& ietaE,
0322                                  std::vector<int>& ietaW,
0323                                  std::vector<int>& iphiN,
0324                                  std::vector<int>& iphiS,
0325                                  std::vector<CaloDirection>& dir,
0326                                  const CaloSubdetectorTopology* barrelTopo,
0327                                  const CaloSubdetectorTopology* endcapTopo,
0328                                  const EcalBarrelGeometry* barrelGeom,
0329                                  const EcalEndcapGeometry* endcapGeom,
0330                                  bool debug,
0331                                  bool ignoreTransition) {
0332     if (debug) {
0333       edm::LogVerbatim("IsoTrack") << "newECALIdNS::Add columns of cells for " << (dets.size() - last)
0334                                    << " cells (last) " << last;
0335       for (unsigned int i1 = last; i1 < dets.size(); i1++) {
0336         edm::LogVerbatim("IsoTrack") << spr::debugEcalDet(i1, dets[i1]).str() << " along " << dir[i1] << " # "
0337                                      << iphiN[i1] << "|" << iphiS[i1];
0338       }
0339     }
0340     std::vector<DetId> vdets;
0341     std::vector<CaloDirection> dirs;
0342     std::vector<int> jetaE, jetaW, jphiN, jphiS;
0343     vdets.insert(vdets.end(), dets.begin(), dets.end());
0344     dirs.insert(dirs.end(), dir.begin(), dir.end());
0345     jetaE.insert(jetaE.end(), ietaE.begin(), ietaE.end());
0346     jetaW.insert(jetaW.end(), ietaW.begin(), ietaW.end());
0347     jphiN.insert(jphiN.end(), iphiN.begin(), iphiN.end());
0348     jphiS.insert(jphiS.end(), iphiS.begin(), iphiS.end());
0349     std::vector<DetId> vdetE, vdetW;
0350     if (last == 0) {
0351       unsigned int ndet = vdets.size();
0352       std::vector<CaloDirection> dirE(ndet, EAST), dirW(ndet, WEST);
0353       vdetE = spr::newECALIdEW(
0354           dets, last, ietaE, ietaW, dirE, barrelTopo, endcapTopo, barrelGeom, endcapGeom, debug, ignoreTransition);
0355       vdetW = spr::newECALIdEW(
0356           dets, last, ietaE, ietaW, dirW, barrelTopo, endcapTopo, barrelGeom, endcapGeom, debug, ignoreTransition);
0357       for (unsigned int i1 = 0; i1 < vdetW.size(); i1++) {
0358         if (std::count(vdets.begin(), vdets.end(), vdetW[i1]) == 0) {
0359           vdets.push_back(vdetW[i1]);
0360           dirs.push_back(dir[0]);
0361           jetaE.push_back(0);
0362           jetaW.push_back(0);
0363           jphiN.push_back(iphiN[0]);
0364           jphiS.push_back(iphiS[0]);
0365         }
0366       }
0367       for (unsigned int i1 = 0; i1 < vdetE.size(); i1++) {
0368         if (std::count(vdets.begin(), vdets.end(), vdetE[i1]) == 0) {
0369           vdets.push_back(vdetE[i1]);
0370           dirs.push_back(dir[0]);
0371           jetaE.push_back(0);
0372           jetaW.push_back(0);
0373           jphiN.push_back(iphiN[0]);
0374           jphiS.push_back(iphiS[0]);
0375         }
0376       }
0377       if (debug) {
0378         edm::LogVerbatim("IsoTrack") << "newECALIdNS::With Added cells along E/W results a set of "
0379                                      << (vdets.size() - dets.size()) << " new  cells";
0380         spr::debugEcalDets(dets.size(), vdets, dirs);
0381       }
0382     }
0383 
0384     unsigned int last0 = vdets.size();
0385     std::vector<DetId> vdetnew;
0386     std::vector<CaloDirection> dirnew;
0387     std::vector<int> kphiN, kphiS, ketaE, ketaW;
0388     int kphi = 0;
0389     for (unsigned int i1 = last; i1 < dets.size(); i1++) {
0390       int iphi = iphiS[i1];
0391       if (dir[i1] == NORTH)
0392         iphi = iphiN[i1];
0393       if (iphi > 0) {
0394         std::vector<DetId> vdetn(1);
0395         std::vector<CaloDirection> dirn(1);
0396         std::vector<CaloDirection> dirnE(1, EAST), dirnW(1, WEST);
0397         int flag = 0;
0398         std::vector<DetId> cells;
0399         spr::simpleMove(
0400             dets[i1], dir[i1], barrelTopo, endcapTopo, barrelGeom, endcapGeom, cells, flag, debug, ignoreTransition);
0401         iphi--;
0402         if (iphi > kphi)
0403           kphi = iphi;
0404         if (dir[i1] == NORTH)
0405           jphiN[i1] = iphi;
0406         else
0407           jphiS[i1] = iphi;
0408         if (flag != 0) {
0409           if (std::count(vdets.begin(), vdets.end(), cells[0]) == 0) {
0410             int kfiN = iphiN[i1];
0411             int kfiS = iphiS[i1];
0412             vdetn[0] = cells[0];
0413             vdetnew.push_back(vdetn[0]);
0414             dirn[0] = dir[i1];
0415             if (dir[i1] == NORTH)
0416               kfiN = iphi;
0417             else
0418               kfiS = iphi;
0419             if (flag < 0) {
0420               int ktmp = kfiS;
0421               kfiS = kfiN;
0422               kfiN = ktmp;
0423               if (dirn[0] == NORTH)
0424                 dirn[0] = SOUTH;
0425               else
0426                 dirn[0] = NORTH;
0427             }
0428             dirnew.push_back(dirn[0]);
0429             kphiN.push_back(kfiN);
0430             ketaE.push_back(ietaE[i1]);
0431             kphiS.push_back(kfiS);
0432             ketaW.push_back(ietaW[i1]);
0433             std::vector<int> ietE(1, ietaE[i1]), ietW(1, ietaW[i1]);
0434             vdetE = spr::newECALIdEW(
0435                 vdetn, 0, ietE, ietW, dirnE, barrelTopo, endcapTopo, barrelGeom, endcapGeom, debug, ignoreTransition);
0436             vdetW = spr::newECALIdEW(
0437                 vdetn, 0, ietE, ietW, dirnW, barrelTopo, endcapTopo, barrelGeom, endcapGeom, debug, ignoreTransition);
0438             for (unsigned int i2 = 0; i2 < vdetW.size(); i2++) {
0439               if (std::count(vdets.begin(), vdets.end(), vdetW[i2]) == 0 &&
0440                   std::count(vdetnew.begin(), vdetnew.end(), vdetW[i2]) == 0) {
0441                 vdets.push_back(vdetW[i2]);
0442                 dirs.push_back(dirn[0]);
0443                 jetaE.push_back(0);
0444                 jphiN.push_back(kfiN);
0445                 jetaW.push_back(0);
0446                 jphiS.push_back(kfiS);
0447               }
0448             }
0449             for (unsigned int i2 = 0; i2 < vdetE.size(); i2++) {
0450               if (std::count(vdets.begin(), vdets.end(), vdetE[i2]) == 0 &&
0451                   std::count(vdetnew.begin(), vdetnew.end(), vdetE[i2]) == 0) {
0452                 vdets.push_back(vdetE[i2]);
0453                 dirs.push_back(dirn[0]);
0454                 jetaE.push_back(0);
0455                 jphiN.push_back(kfiN);
0456                 jetaW.push_back(0);
0457                 jphiS.push_back(kfiS);
0458               }
0459             }
0460           }
0461         }
0462       }
0463     }
0464     last = vdets.size();
0465     for (unsigned int i2 = 0; i2 < vdetnew.size(); i2++) {
0466       if (std::count(vdets.begin(), vdets.end(), vdetnew[i2]) == 0) {
0467         vdets.push_back(vdetnew[i2]);
0468         dirs.push_back(dirnew[i2]);
0469         jetaE.push_back(ketaE[i2]);
0470         jetaW.push_back(ketaW[i2]);
0471         jphiN.push_back(kphiN[i2]);
0472         jphiS.push_back(kphiS[i2]);
0473       }
0474     }
0475     if (debug) {
0476       edm::LogVerbatim("IsoTrack") << "newECALIdNS::Addition results a set of " << (vdets.size() - last0)
0477                                    << " new  cells (last " << last0 << ", iphi " << kphi << ")";
0478       for (unsigned int i1 = last0; i1 < vdets.size(); i1++) {
0479         edm::LogVerbatim("IsoTrack") << spr::debugEcalDet(i1, vdets[i1]).str() << " along " << dirs[i1] << " iphi "
0480                                      << jphiN[i1] << "|" << jphiS[i1] << " ieta " << jetaE[i1] << "|" << jetaW[i1];
0481       }
0482     }
0483     last0 = last;
0484 
0485     if (kphi > 0) {
0486       last = last0;
0487       return spr::newECALIdNS(vdets,
0488                               last,
0489                               jetaE,
0490                               jetaW,
0491                               jphiN,
0492                               jphiS,
0493                               dirs,
0494                               barrelTopo,
0495                               endcapTopo,
0496                               barrelGeom,
0497                               endcapGeom,
0498                               debug,
0499                               ignoreTransition);
0500     } else {
0501       if (debug) {
0502         edm::LogVerbatim("IsoTrack") << "newECALIdNS::Final list consists of " << vdets.size() << " cells";
0503         spr::debugEcalDets(0, vdets);
0504       }
0505       return vdets;
0506     }
0507   }
0508 
0509   std::vector<DetId> newECALIdEW(std::vector<DetId>& dets,
0510                                  unsigned int last,
0511                                  int ieta,
0512                                  std::vector<CaloDirection>& dir,
0513                                  const CaloSubdetectorTopology* barrelTopo,
0514                                  const CaloSubdetectorTopology* endcapTopo,
0515                                  const EcalBarrelGeometry* barrelGeom,
0516                                  const EcalEndcapGeometry* endcapGeom,
0517                                  bool debug,
0518                                  bool ignoreTransition) {
0519     if (debug) {
0520       edm::LogVerbatim("IsoTrack") << "newECALIdEW::Add " << ieta << " rows of cells for " << last << ":" << dets.size()
0521                                    << ":" << (dets.size() - last) << " cells" << std::endl;
0522       spr::debugEcalDets(last, dets, dir);
0523     }
0524     std::vector<DetId> vdets;
0525     vdets.clear();
0526     std::vector<CaloDirection> dirs;
0527     dirs.clear();
0528     vdets.insert(vdets.end(), dets.begin(), dets.end());
0529     dirs.insert(dirs.end(), dir.begin(), dir.end());
0530 
0531     if (ieta > 0) {
0532       for (unsigned int i1 = last; i1 < dets.size(); i1++) {
0533         int flag = 0;
0534         std::vector<DetId> cells;
0535         spr::simpleMove(
0536             dets[i1], dir[i1], barrelTopo, endcapTopo, barrelGeom, endcapGeom, cells, flag, debug, ignoreTransition);
0537         if (flag != 0) {
0538           if (std::count(vdets.begin(), vdets.end(), cells[0]) == 0) {
0539             CaloDirection dirn = dir[i1];
0540             if (flag < 0) {
0541               if (dirn == EAST)
0542                 dirn = WEST;
0543               else
0544                 dirn = EAST;
0545             }
0546             vdets.push_back(cells[0]);
0547             dirs.push_back(dirn);
0548           }
0549         }
0550       }
0551       ieta--;
0552     }
0553 
0554     if (debug) {
0555       edm::LogVerbatim("IsoTrack") << "newECALIdEW::Addition results a set of " << (vdets.size() - dets.size())
0556                                    << " new  cells";
0557       spr::debugEcalDets(dets.size(), vdets, dirs);
0558     }
0559     if (ieta > 0) {
0560       last = dets.size();
0561       return spr::newECALIdEW(
0562           vdets, last, ieta, dirs, barrelTopo, endcapTopo, barrelGeom, endcapGeom, debug, ignoreTransition);
0563     } else {
0564       if (debug) {
0565         edm::LogVerbatim("IsoTrack") << "newECALIdEW::Final list (EW) consists of " << vdets.size() << " cells";
0566         spr::debugEcalDets(0, vdets);
0567       }
0568       return vdets;
0569     }
0570   }
0571 
0572   std::vector<DetId> newECALIdEW(std::vector<DetId>& dets,
0573                                  unsigned int last,
0574                                  std::vector<int>& ietaE,
0575                                  std::vector<int>& ietaW,
0576                                  std::vector<CaloDirection>& dir,
0577                                  const CaloSubdetectorTopology* barrelTopo,
0578                                  const CaloSubdetectorTopology* endcapTopo,
0579                                  const EcalBarrelGeometry* barrelGeom,
0580                                  const EcalEndcapGeometry* endcapGeom,
0581                                  bool debug,
0582                                  bool ignoreTransition) {
0583     if (debug) {
0584       edm::LogVerbatim("IsoTrack") << "newECALIdEW::Add " << ietaE[0] << "|" << ietaW[0] << " rows of cells for "
0585                                    << (dets.size() - last) << " cells (last " << last << ")";
0586       spr::debugEcalDets(last, dets, dir);
0587     }
0588     std::vector<DetId> vdets;
0589     vdets.insert(vdets.end(), dets.begin(), dets.end());
0590     std::vector<CaloDirection> dirs;
0591     dirs.insert(dirs.end(), dir.begin(), dir.end());
0592     std::vector<int> jetaE, jetaW;
0593     jetaE.insert(jetaE.end(), ietaE.begin(), ietaE.end());
0594     jetaW.insert(jetaW.end(), ietaW.begin(), ietaW.end());
0595     int keta = 0;
0596     for (unsigned int i1 = last; i1 < dets.size(); i1++) {
0597       int ieta = ietaW[i1];
0598       if (dir[i1] == EAST)
0599         ieta = ietaE[i1];
0600       if (ieta > 0) {
0601         int flag = 0;
0602         std::vector<DetId> cells;
0603         spr::simpleMove(
0604             dets[i1], dir[i1], barrelTopo, endcapTopo, barrelGeom, endcapGeom, cells, flag, debug, ignoreTransition);
0605         ieta--;
0606         if (ieta > keta)
0607           keta = ieta;
0608         if (dir[i1] == EAST)
0609           jetaE[i1] = ieta;
0610         else
0611           jetaW[i1] = ieta;
0612         if (flag != 0) {
0613           if (std::count(vdets.begin(), vdets.end(), cells[0]) == 0) {
0614             vdets.push_back(cells[0]);
0615             CaloDirection dirn = dir[i1];
0616             int ketaE = ietaE[i1];
0617             int ketaW = ietaW[i1];
0618             if (dirn == EAST)
0619               ketaE = ieta;
0620             else
0621               ketaW = ieta;
0622             if (flag < 0) {
0623               int ktmp = ketaW;
0624               ketaW = ketaE;
0625               ketaE = ktmp;
0626               if (dirn == EAST)
0627                 dirn = WEST;
0628               else
0629                 dirn = EAST;
0630             }
0631             dirs.push_back(dirn);
0632             jetaE.push_back(ketaE);
0633             jetaW.push_back(ketaW);
0634           }
0635         }
0636       }
0637     }
0638 
0639     if (debug) {
0640       edm::LogVerbatim("IsoTrack") << "newECALIdEW::Addition results a set of " << (vdets.size() - dets.size())
0641                                    << " new  cells (last " << dets.size() << ", ieta " << keta << ")";
0642       spr::debugEcalDets(dets.size(), vdets);
0643     }
0644     if (keta > 0) {
0645       last = dets.size();
0646       return spr::newECALIdEW(
0647           vdets, last, jetaE, jetaW, dirs, barrelTopo, endcapTopo, barrelGeom, endcapGeom, debug, ignoreTransition);
0648     } else {
0649       if (debug) {
0650         edm::LogVerbatim("IsoTrack") << "newECALIdEW::Final list (EW) consists of " << vdets.size() << " cells";
0651         spr::debugEcalDets(0, vdets);
0652       }
0653       return vdets;
0654     }
0655   }
0656 
0657   void simpleMove(DetId& det,
0658                   const CaloDirection& dir,
0659                   const CaloSubdetectorTopology* barrelTopo,
0660                   const CaloSubdetectorTopology* endcapTopo,
0661                   const EcalBarrelGeometry* barrelGeom,
0662                   const EcalEndcapGeometry* endcapGeom,
0663                   std::vector<DetId>& cells,
0664                   int& ok,
0665                   bool debug,
0666                   bool ignoreTransition) {
0667     DetId cell;
0668     ok = 0;
0669     if (det.subdetId() == EcalBarrel) {
0670       EBDetId detId = det;
0671       std::vector<DetId> neighbours = barrelTopo->getNeighbours(detId, dir);
0672       if (!neighbours.empty() && !neighbours[0].null()) {
0673         cells.push_back(neighbours[0]);
0674         cell = neighbours[0];
0675         ok = 1;
0676       } else {
0677         const int ietaAbs(detId.ietaAbs());  // abs value of ieta
0678         if (EBDetId::MAX_IETA == ietaAbs && (!ignoreTransition) && endcapGeom) {
0679           // get ee nbrs for for end of barrel crystals
0680           const EcalBarrelGeometry::OrderedListOfEEDetId& ol(*barrelGeom->getClosestEndcapCells(detId));
0681           // take closest neighbour on the other side, that is in the endcap
0682           cell = *(ol.begin());
0683           neighbours = endcapTopo->getNeighbours(cell, dir);
0684           if (!neighbours.empty() && !neighbours[0].null())
0685             ok = 1;
0686           else
0687             ok = -1;
0688           for (EcalBarrelGeometry::OrderedListOfEEDetId::const_iterator iptr = ol.begin(); iptr != ol.end(); ++iptr)
0689             cells.push_back(*iptr);
0690         }
0691       }
0692     } else if (det.subdetId() == EcalEndcap && endcapGeom) {
0693       EEDetId detId = det;
0694       std::vector<DetId> neighbours = endcapTopo->getNeighbours(detId, dir);
0695       if (!neighbours.empty() && !neighbours[0].null()) {
0696         cells.push_back(neighbours[0]);
0697         cell = neighbours[0];
0698         ok = 1;
0699       } else {
0700         // are we on the outer ring ?
0701         const int iphi(detId.iPhiOuterRing());
0702         //  int isc = detId.isc();
0703         if (iphi != 0 && (!ignoreTransition)) {
0704           // get eb nbrs for for end of endcap crystals
0705           const EcalEndcapGeometry::OrderedListOfEBDetId& ol(*endcapGeom->getClosestBarrelCells(detId));
0706           // take closest neighbour on the other side, that is in the barrel.
0707           cell = *(ol.begin());
0708           neighbours = barrelTopo->getNeighbours(cell, dir);
0709           if (!neighbours.empty() && !neighbours[0].null())
0710             ok = 1;
0711           else
0712             ok = -1;
0713           for (EcalEndcapGeometry::OrderedListOfEBDetId::const_iterator iptr = ol.begin(); iptr != ol.end(); ++iptr)
0714             cells.push_back(*iptr);
0715         }
0716       }
0717     }
0718     if (debug) {
0719       std::ostringstream st1;
0720       for (unsigned int i1 = 0; i1 < cells.size(); ++i1)
0721         st1 << " " << std::hex << cells[0]() << std::dec;
0722       edm::LogVerbatim("IsoTrack") << "simpleMove:: Move DetId 0x" << std::hex << det() << std::dec << " along " << dir
0723                                    << " to get 0x" << std::hex << cell() << std::dec << " with flag " << ok << " # "
0724                                    << cells.size() << st1.str();
0725     }
0726   }
0727 
0728   void extraIds(const DetId& det,
0729                 std::vector<DetId>& dets,
0730                 int ietaE,
0731                 int ietaW,
0732                 int iphiN,
0733                 int iphiS,
0734                 const EcalBarrelGeometry* barrelGeom,
0735                 const EcalEndcapGeometry* endcapGeom,
0736                 std::vector<DetId>& cells,
0737                 bool debug) {
0738     if (det.subdetId() == EcalBarrel) {
0739       EBDetId id = det;
0740       if (debug)
0741         edm::LogVerbatim("IsoTrack") << "extraIds::Cell " << id << " rows " << ietaW << "|" << ietaE << " columns "
0742                                      << iphiS << "|" << iphiN;
0743       int etaC = id.ietaAbs();
0744       int phiC = id.iphi();
0745       int zsid = id.zside();
0746       for (int eta = -ietaW; eta <= ietaE; ++eta) {
0747         for (int phi = -iphiS; phi <= iphiN; ++phi) {
0748           int iphi = phiC + phi;
0749           if (iphi < 0)
0750             iphi += 360;
0751           else if (iphi > 360)
0752             iphi -= 360;
0753           int ieta = zsid * (etaC + eta);
0754           if (EBDetId::validDetId(ieta, iphi)) {
0755             id = EBDetId(ieta, iphi);
0756             if (barrelGeom->present(id)) {
0757               if (std::count(dets.begin(), dets.end(), (DetId)id) == 0) {
0758                 cells.push_back((DetId)id);
0759               }
0760             }
0761           }
0762         }
0763       }
0764     } else if (det.subdetId() == EcalEndcap && endcapGeom) {
0765       EEDetId id = det;
0766       if (debug)
0767         edm::LogVerbatim("IsoTrack") << "extraIds::Cell " << id << " rows " << ietaW << "|" << ietaE << " columns "
0768                                      << iphiS << "|" << iphiN;
0769       int ixC = id.ix();
0770       int iyC = id.iy();
0771       int zsid = id.zside();
0772       for (int kx = -ietaW; kx <= ietaE; ++kx) {
0773         for (int ky = -iphiS; ky <= iphiN; ++ky) {
0774           int ix = ixC + kx;
0775           int iy = iyC + ky;
0776           if (EEDetId::validDetId(ix, iy, zsid)) {
0777             id = EEDetId(ix, iy, zsid);
0778             if (endcapGeom->present(id)) {
0779               if (std::count(dets.begin(), dets.end(), (DetId)id) == 0) {
0780                 cells.push_back((DetId)id);
0781               }
0782             }
0783           }
0784         }
0785       }
0786     }
0787 
0788     if (debug) {
0789       edm::LogVerbatim("IsoTrack") << "extraIds:: finds " << cells.size() << " new cells";
0790       spr::debugEcalDets(0, cells);
0791     }
0792   }
0793 }  // namespace spr
[ Source navigation ] [ Diff markup ] [ Identifier search ] [ general search ]

This page was automatically generated by the 2.2.1 LXR engine.

The LXR team
Valid CSS 2.1! Valid HTML 4.01!