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 }
1.7.3