CMS 3D CMS Logo

/data/refman/pasoursint/CMSSW_5_3_3/src/Calibration/IsolatedParticles/src/MatrixHCALDetIds.cc

Go to the documentation of this file.
00001 #include "DataFormats/HcalDetId/interface/HcalDetId.h"
00002 
00003 #include "Calibration/IsolatedParticles/interface/MatrixHCALDetIds.h"
00004 #include "Calibration/IsolatedParticles/interface/FindDistCone.h"
00005 #include "Calibration/IsolatedParticles/interface/DebugInfo.h"
00006 
00007 #include<algorithm>
00008 #include<iostream>
00009 
00010 namespace spr{
00011 
00012   std::vector<DetId> matrixHCALIds(std::vector<DetId>& dets, 
00013                                    const HcalTopology* topology, int ieta, 
00014                                    int iphi, bool includeHO, bool debug) {
00015 
00016     if (debug) {
00017       std::cout << "matrixHCALIds::Add " << ieta << " rows and " << iphi 
00018                 << " columns of cells for " << dets.size() << " cells" 
00019                 << std::endl;
00020       spr::debugHcalDets(0, dets);
00021     }
00022 
00023     std::vector<DetId> vdetN = spr::newHCALIdNS(dets, 0, topology, true,  ieta,
00024                                                iphi, debug);
00025     std::vector<DetId> vdetS = spr::newHCALIdNS(dets, 0, topology, false, ieta,
00026                                                 iphi, debug);
00027     for (unsigned int i1=0; i1<vdetS.size(); i1++) {
00028       if (std::count(vdetN.begin(),vdetN.end(),vdetS[i1]) == 0)
00029         vdetN.push_back(vdetS[i1]);
00030     }
00031 
00032     vdetS = spr::matrixHCALIdsDepth(vdetN, topology, includeHO, debug);
00033 
00034     if (debug) {
00035       std::cout << "matrixHCALIds::Total number of cells found is " 
00036                 << vdetS.size() << std::endl;
00037       spr::debugHcalDets(0, vdetS);
00038     }
00039     return vdetS;
00040   }
00041 
00042   std::vector<DetId> matrixHCALIds(const DetId& det, const CaloGeometry* geo,
00043                                    const HcalTopology* topology, double dR, 
00044                                    const GlobalVector& trackMom, bool includeHO,
00045                                    bool debug) {
00046  
00047     HcalDetId   hcdet = HcalDetId(det);
00048     GlobalPoint core  = geo->getPosition(hcdet);
00049     std::vector<DetId> dets, vdetx;
00050     dets.push_back(det);
00051     int ietaphi = (int)(dR/15.0)+1;
00052     std::vector<DetId> vdets = spr::matrixHCALIds(dets, topology, ietaphi, 
00053                                                   ietaphi, includeHO, debug);
00054     for (unsigned int i=0; i<vdets.size(); ++i) {
00055       HcalDetId   hcdet  = HcalDetId(vdets[i]);
00056       GlobalPoint rpoint = geo->getPosition(hcdet);
00057       if (spr::getDistInPlaneTrackDir(core, trackMom, rpoint) < dR) {
00058         vdetx.push_back(vdets[i]);
00059       }
00060     }
00061 
00062     if (debug) {
00063       std::cout << "matrixHCALIds::Final List of cells for dR " << dR
00064                 << " is with " << vdetx.size() << " from original list of "
00065                 << vdets.size() << " cells" << std::endl;
00066       spr::debugHcalDets(0, vdetx);
00067     }
00068     return vdetx;
00069  }
00070 
00071   std::vector<DetId> matrixHCALIds(std::vector<DetId>& dets, 
00072                                    const HcalTopology* topology, int ietaE, 
00073                                    int ietaW,int iphiN,int iphiS, 
00074                                    bool includeHO, bool debug) {
00075 
00076     if (debug) {
00077       std::cout << "matrixHCALIds::Add " <<ietaE << "|" <<ietaW << " rows and "
00078                 << iphiN << "|" << iphiS << " columns of cells for " 
00079                 << dets.size() << " cells" << std::endl;
00080       spr::debugHcalDets(0, dets);
00081     }
00082 
00083     std::vector<DetId> vdetN = spr::newHCALIdNS(dets, 0, topology, true, ietaE,
00084                                                 ietaW, iphiN, iphiS, debug);
00085     std::vector<DetId> vdetS = spr::newHCALIdNS(dets, 0, topology, false,ietaE,
00086                                                 ietaW, iphiN, iphiS, debug);
00087     for (unsigned int i1=0; i1<vdetS.size(); i1++) {
00088       if (std::count(vdetN.begin(),vdetN.end(),vdetS[i1]) == 0)
00089         vdetN.push_back(vdetS[i1]);
00090     }
00091 
00092     vdetS = spr::matrixHCALIdsDepth(vdetN, topology, includeHO, debug);
00093 
00094     if (debug) {
00095       std::cout << "matrixHCALIds::Total number of cells found is " 
00096                 << vdetS.size() << std::endl;
00097       spr::debugHcalDets(0, vdetS);
00098     }
00099     return vdetS;
00100   }
00101 
00102   std::vector<DetId> newHCALIdNS(std::vector<DetId>& dets, unsigned int last,
00103                                  const HcalTopology* topology, bool shiftNorth,
00104                                  int ieta, int iphi, bool debug) {
00105 
00106     if (debug) {
00107       std::cout << "newHCALIdNS::Add " << iphi << " columns of cells along " 
00108                 << shiftNorth << " for " << (dets.size()-last) << " cells" 
00109                 << std::endl;
00110       spr::debugHcalDets(last, dets);
00111     }
00112 
00113     std::vector<DetId> vdets;
00114     vdets.insert(vdets.end(), dets.begin(), dets.end());
00115     std::vector<DetId> vdetE, vdetW;
00116     if (last == 0) {
00117       vdetE = spr::newHCALIdEW(dets, last, topology, true,  ieta, debug);
00118       vdetW = spr::newHCALIdEW(dets, last, topology, false, ieta, debug);
00119       for (unsigned int i1=0; i1<vdetW.size(); i1++) {
00120         if (std::count(vdets.begin(),vdets.end(),vdetW[i1]) == 0)
00121           vdets.push_back(vdetW[i1]);
00122       }
00123       for (unsigned int i1=0; i1<vdetE.size(); i1++) {
00124         if (std::count(vdets.begin(),vdets.end(),vdetE[i1]) == 0)
00125           vdets.push_back(vdetE[i1]);
00126       }
00127       if (debug) {
00128         std::cout <<"newHCALIdNS::With Added cells along E/W results a set of "
00129                   << (vdets.size()-dets.size()) << " new  cells" << std::endl;
00130         spr::debugHcalDets(dets.size(), vdets);
00131       }
00132     }
00133     unsigned int last0 = vdets.size();
00134     if (iphi > 0) {
00135       std::vector<DetId> vdetnew;
00136       for (unsigned int i1=last; i1<dets.size(); i1++) {
00137         std::vector<DetId> vdet;
00138         if (shiftNorth) vdet = topology->north(dets[i1]);
00139         else            vdet = topology->south(dets[i1]);
00140         for (unsigned int i2=0; i2<vdet.size(); i2++) {
00141           if (std::count(vdets.begin(),vdets.end(),vdet[i2]) == 0)
00142             vdetnew.push_back(vdet[i2]);
00143         }
00144       }
00145       iphi--;
00146       vdetE = spr::newHCALIdEW(vdetnew, 0, topology, true,  ieta, debug);
00147       vdetW = spr::newHCALIdEW(vdetnew, 0, topology, false, ieta, debug);
00148       for (unsigned int i2=0; i2<vdetW.size(); i2++) {
00149         if (std::count(vdets.begin(),vdets.end(),vdetW[i2]) == 0 &&
00150             std::count(vdetnew.begin(),vdetnew.end(),vdetW[i2]) == 0)
00151           vdets.push_back(vdetW[i2]);
00152       }
00153       for (unsigned int i2=0; i2<vdetE.size(); i2++) {
00154         if (std::count(vdets.begin(),vdets.end(),vdetE[i2]) == 0 &&
00155             std::count(vdetnew.begin(),vdetnew.end(),vdetE[i2]) == 0)
00156           vdets.push_back(vdetE[i2]);
00157       }
00158       last = vdets.size();
00159       vdets.insert(vdets.end(), vdetnew.begin(), vdetnew.end());
00160       if (debug) {
00161         std::cout << "newHCALIdNS::Addition results a set of " 
00162                   << (vdets.size()-last0)  << " new  cells" << std::endl;
00163         spr::debugHcalDets(last0, vdets);
00164       }
00165       last0 = last;
00166     }
00167 
00168     if (iphi > 0) {
00169       last = last0;
00170       return spr::newHCALIdNS(vdets,last,topology,shiftNorth,ieta,iphi,debug);
00171     } else {
00172       if (debug) {
00173         std::cout << "newHCALIdNS::Final list consists of " << vdets.size()
00174                   << " cells" << std::endl;
00175         spr::debugHcalDets(0, vdets);
00176       }
00177       return vdets;
00178     }
00179   }
00180 
00181   std::vector<DetId> newHCALIdNS(std::vector<DetId>& dets, unsigned int last,
00182                                  const HcalTopology* topology, bool shiftNorth,
00183                                  int ietaE, int ietaW, int iphiN, int iphiS,
00184                                  bool debug) {
00185 
00186     if (debug) {
00187       std::cout << "newHCALIdNS::Add " << iphiN << "|" << iphiS
00188                 << " columns of cells along " << shiftNorth << " for " 
00189                 << (dets.size()-last) << " cells" << std::endl;
00190       spr::debugHcalDets(last, dets);
00191     }
00192 
00193     std::vector<DetId> vdets;
00194     vdets.insert(vdets.end(), dets.begin(), dets.end());
00195     std::vector<DetId> vdetE, vdetW;
00196     if (last == 0) {
00197       vdetE = spr::newHCALIdEW(dets,last, topology, true,  ietaE,ietaW, debug);
00198       vdetW = spr::newHCALIdEW(dets,last, topology, false, ietaE,ietaW, debug);
00199       for (unsigned int i1=0; i1<vdetW.size(); i1++) {
00200         if (std::count(vdets.begin(),vdets.end(),vdetW[i1]) == 0)
00201           vdets.push_back(vdetW[i1]);
00202       }
00203       for (unsigned int i1=0; i1<vdetE.size(); i1++) {
00204         if (std::count(vdets.begin(),vdets.end(),vdetE[i1]) == 0)
00205           vdets.push_back(vdetE[i1]);
00206       }
00207       if (debug) {
00208         std::cout <<"newHCALIdNS::With Added cells along E/W results a set of "
00209                   << (vdets.size()-dets.size()) << " new  cells" << std::endl;
00210         spr::debugHcalDets(dets.size(), vdets);
00211       }
00212     }
00213     unsigned int last0 = vdets.size();
00214     int iphi = iphiS;
00215     if (shiftNorth) iphi = iphiN;
00216     if (iphi > 0) {
00217       std::vector<DetId> vdetnew;
00218       for (unsigned int i1=last; i1<dets.size(); i1++) {
00219         std::vector<DetId> vdet;
00220         if (shiftNorth) vdet = topology->north(dets[i1]);
00221         else            vdet = topology->south(dets[i1]);
00222         for (unsigned int i2=0; i2<vdet.size(); i2++) {
00223           if (std::count(vdets.begin(),vdets.end(),vdet[i2]) == 0)
00224             vdetnew.push_back(vdet[i2]);
00225         }
00226       }
00227       iphi--;
00228       vdetE = spr::newHCALIdEW(vdetnew,0, topology, true,  ietaE,ietaW, debug);
00229       vdetW = spr::newHCALIdEW(vdetnew,0, topology, false, ietaE,ietaW, debug);
00230       for (unsigned int i2=0; i2<vdetW.size(); i2++) {
00231         if (std::count(vdets.begin(),vdets.end(),vdetW[i2]) == 0 &&
00232             std::count(vdetnew.begin(),vdetnew.end(),vdetW[i2]) == 0)
00233           vdets.push_back(vdetW[i2]);
00234       }
00235       for (unsigned int i2=0; i2<vdetE.size(); i2++) {
00236         if (std::count(vdets.begin(),vdets.end(),vdetE[i2]) == 0 &&
00237             std::count(vdetnew.begin(),vdetnew.end(),vdetE[i2]) == 0)
00238           vdets.push_back(vdetE[i2]);
00239       }
00240       last = vdets.size();
00241       vdets.insert(vdets.end(), vdetnew.begin(), vdetnew.end());
00242       if (debug) {
00243         std::cout << "newHCALIdNS::Addition results a set of " 
00244                   << (vdets.size()-last0)  << " new  cells" << std::endl;
00245         spr::debugHcalDets(last0, vdets);
00246       }
00247       last0 = last;
00248     }
00249     if (shiftNorth) iphiN = iphi;
00250     else            iphiS = iphi;
00251 
00252     if (iphi > 0) {
00253       last = last0;
00254       return spr::newHCALIdNS(vdets,last,topology,shiftNorth,ietaE,ietaW,
00255                               iphiN,iphiS,debug);
00256     } else {
00257       if (debug) {
00258         std::cout << "newHCALIdNS::Final list consists of " << vdets.size()
00259                   << " cells" << std::endl;
00260         spr::debugHcalDets(0, vdets);
00261       }
00262       return vdets;
00263     }
00264   }
00265 
00266   std::vector<DetId> newHCALIdEW(std::vector<DetId>& dets, unsigned int last,
00267                                  const HcalTopology* topology, bool shiftEast,
00268                                  int ieta, bool debug) {
00269 
00270     if (debug) {
00271       std::cout << "newHCALIdEW::Add " << ieta << " rows of cells along " 
00272                 << shiftEast << " for " << (dets.size()-last) << " cells" 
00273                 << std::endl;
00274       spr::debugHcalDets(last, dets);
00275     }
00276 
00277     std::vector<DetId> vdets;
00278     vdets.insert(vdets.end(), dets.begin(), dets.end());
00279     if (ieta > 0) {
00280       for (unsigned int i1=last; i1<dets.size(); i1++) {
00281         std::vector<DetId> vdet;
00282         if (shiftEast) vdet = topology->east(dets[i1]);
00283         else           vdet = topology->west(dets[i1]);
00284         for (unsigned int i2=0; i2<vdet.size(); i2++) {
00285           if (std::count(vdets.begin(),vdets.end(),vdet[i2]) == 0)
00286             vdets.push_back(vdet[i2]);
00287         }
00288       }
00289       ieta--;
00290     }
00291     
00292     if (debug) {
00293       std::cout << "newHCALIdEW::Addition results a set of " 
00294                 << (vdets.size()-dets.size()) << " new  cells" << std::endl;
00295       spr::debugHcalDets(dets.size(), vdets);
00296     }
00297 
00298     if (ieta > 0) {
00299       last = dets.size();
00300       return spr::newHCALIdEW(vdets, last, topology, shiftEast, ieta, debug);
00301     } else {
00302       if (debug) {
00303         std::cout << "newHCALIdEW::Final list (EW) consists of " <<vdets.size()
00304                   << " cells" << std::endl;
00305         spr::debugHcalDets(0, vdets);
00306       }
00307       return vdets;
00308     }
00309   }
00310 
00311   std::vector<DetId> newHCALIdEW(std::vector<DetId>& dets, unsigned int last,
00312                                  const HcalTopology* topology, bool shiftEast,
00313                                  int ietaE, int ietaW, bool debug) {
00314 
00315     if (debug) {
00316       std::cout << "newHCALIdEW::Add " << ietaE << "|" << ietaW
00317                 << " rows of cells along " << shiftEast << " for " 
00318                 << (dets.size()-last) << " cells" << std::endl;
00319       spr::debugHcalDets(last, dets);
00320     }
00321 
00322     int ieta = ietaW;
00323     if (shiftEast) ieta = ietaE;
00324     std::vector<DetId> vdets;
00325     vdets.insert(vdets.end(), dets.begin(), dets.end());
00326     if (ieta > 0) {
00327       for (unsigned int i1=last; i1<dets.size(); i1++) {
00328         std::vector<DetId> vdet;
00329         if (shiftEast) vdet = topology->east(dets[i1]);
00330         else           vdet = topology->west(dets[i1]);
00331         for (unsigned int i2=0; i2<vdet.size(); i2++) {
00332           if (std::count(vdets.begin(),vdets.end(),vdet[i2]) == 0)
00333             vdets.push_back(vdet[i2]);
00334         }
00335       }
00336       ieta--;
00337     }
00338     if (shiftEast) ietaE = ieta;
00339     else           ietaW = ieta;
00340     
00341     if (debug) {
00342       std::cout << "newHCALIdEW::Addition results a set of " 
00343                 << (vdets.size()-dets.size()) << " new  cells" << std::endl;
00344       spr::debugHcalDets(dets.size(), vdets);
00345     }
00346 
00347     if (ieta > 0) {
00348       last = dets.size();
00349       return spr::newHCALIdEW(vdets,last,topology,shiftEast,ietaE,ietaW,debug);
00350     } else {
00351       if (debug) {
00352         std::cout << "newHCALIdEW::Final list (EW) consists of " <<vdets.size()
00353                   << " cells" << std::endl;
00354         spr::debugHcalDets(0, vdets);
00355       }
00356       return vdets;
00357     }
00358   }
00359 
00360   std::vector<DetId> matrixHCALIdsDepth(std::vector<DetId>& dets, 
00361                                         const HcalTopology* topology, 
00362                                         bool includeHO, bool debug) {
00363 
00364     if (debug) {
00365       std::cout << "matrixHCALIdsDepth::Add cells with higher depths with HO" 
00366                 << "Flag set to " << includeHO << " to existing "
00367                 << dets.size() << " cells" << std::endl;
00368       spr::debugHcalDets(0, dets);
00369     }
00370  
00371     std::vector<DetId> vdets(dets);
00372     for (unsigned int i1=0; i1<dets.size(); i1++) {
00373       HcalDetId vdet = dets[i1];
00374       for (int idepth = 0; idepth < 3; idepth++) {
00375         std::vector<DetId> vUpDetId = topology->up(vdet);
00376         if (vUpDetId.size() != 0) {
00377           if (includeHO || vUpDetId[0].subdetId() != (int)(HcalOuter)) {
00378             int n = std::count(vdets.begin(),vdets.end(),vUpDetId[0]);
00379             if (n == 0) {
00380               if (debug) std::cout << "matrixHCALIdsDepth:: Depth " << idepth 
00381                                    << " " << vdet << " " 
00382                                    << (HcalDetId)vUpDetId[0] << std::endl;
00383               vdets.push_back(vUpDetId[0]);
00384             }
00385           }
00386           vdet = vUpDetId[0];
00387         }
00388       }
00389     }
00390 
00391     if (debug) {
00392       std::cout << "matrixHCALIdsDepth::Final list contains " << vdets.size() 
00393                 << " cells" << std::endl;
00394       spr::debugHcalDets(0, vdets);
00395     }
00396     return vdets;
00397   }
00398 
00399 }