00001 #include "DataFormats/EcalDetId/interface/EBDetId.h"
00002 #include "DataFormats/EcalDetId/interface/EcalSubdetector.h"
00003 #include "DataFormats/EcalDetId/interface/EEDetId.h"
00004
00005 #include "Calibration/IsolatedParticles/interface/MatrixECALDetIds.h"
00006 #include "Calibration/IsolatedParticles/interface/FindDistCone.h"
00007
00008 #include <algorithm>
00009 #include <iostream>
00010
00011 namespace spr{
00012
00013 void matrixECALIds(const DetId& det, int ieta, int iphi,
00014 const CaloGeometry* geo, const CaloTopology* caloTopology,
00015 std::vector<DetId>& vdets, bool debug) {
00016
00017 const CaloSubdetectorTopology *barrelTopo = (caloTopology->getSubdetectorTopology(DetId::Ecal,EcalBarrel));
00018 const CaloSubdetectorTopology *endcapTopo = (caloTopology->getSubdetectorTopology(DetId::Ecal,EcalEndcap));
00019 const EcalBarrelGeometry *barrelGeom = (dynamic_cast< const EcalBarrelGeometry *> (geo->getSubdetectorGeometry(DetId::Ecal,EcalBarrel)));
00020 const EcalEndcapGeometry *endcapGeom = (dynamic_cast< const EcalEndcapGeometry *> (geo->getSubdetectorGeometry(DetId::Ecal,EcalEndcap)));
00021
00022 if (debug) {
00023 std::cout << "matrixECALIds::Add " << ieta << " rows and " << iphi
00024 << " columns of cells for 1 cell" << std::endl;
00025 if (det.subdetId() == EcalBarrel) {
00026 EBDetId id = det;
00027 std::cout << "matrixECALIds::Cell 0x" << std::hex << det() << std::dec
00028 << " " << id << std::endl;
00029 } else if (det.subdetId() == EcalEndcap) {
00030 EEDetId id = det;
00031 std::cout << "matrixECALIds::Cell 0x" << std::hex << det() << std::dec
00032 << " " << id << std::endl;
00033 } else {
00034 std::cout << "matrixECALIds::Cell 0x" << std::hex << det() << std::dec
00035 << " Unknown Type" << std::endl;
00036 }
00037 }
00038
00039 std::vector<DetId> dets(1,det);
00040 std::vector<CaloDirection> dirs(1,NORTH);
00041 vdets = spr::newECALIdNS(dets, 0, ieta,iphi, dirs, *barrelTopo,*endcapTopo,
00042 *barrelGeom,*endcapGeom,debug);
00043 dirs[0] = SOUTH;
00044 std::vector<DetId> vdetS = spr::newECALIdNS(dets, 0, ieta, iphi, dirs,
00045 *barrelTopo,*endcapTopo,
00046 *barrelGeom,*endcapGeom,debug);
00047 for (unsigned int i1=0; i1<vdetS.size(); i1++) {
00048 if (std::count(vdets.begin(),vdets.end(),vdetS[i1]) == 0)
00049 vdets.push_back(vdetS[i1]);
00050 }
00051 unsigned int ndet = (2*ieta+1)*(2*iphi+1);
00052 if (vdets.size() != ndet) {
00053 std::vector<DetId> vdetExtra;
00054 spr::extraIds(det, vdets, ieta, ieta, iphi, iphi,
00055 *barrelGeom, *endcapGeom, vdetExtra, debug);
00056 if (vdetExtra.size() > 0)
00057 vdets.insert(vdets.end(), vdetExtra.begin(), vdetExtra.end());
00058 }
00059
00060 if (debug) {
00061 std::cout << "matrixECALIds::Total number of cells found is "
00062 << vdets.size() << std::endl;
00063 for (unsigned int i1=0; i1<vdets.size(); i1++) {
00064 if (vdets[i1].subdetId() == EcalBarrel) {
00065 EBDetId id = vdets[i1];
00066 std::cout << "matrixECALIds::Cell " << i1 << " 0x" << std::hex
00067 << vdets[i1]() << std::dec << " " << id << std::endl;
00068 } else if (vdets[i1].subdetId() == EcalEndcap) {
00069 EEDetId id = vdets[i1];
00070 std::cout << "matrixECALIds::Cell " << i1 << " 0x" << std::hex
00071 << vdets[i1]() << std::dec << " " << id << std::endl;
00072 } else {
00073 std::cout << "matrixECALIds::Cell " << i1 << " 0x" << std::hex
00074 << vdets[i1]() << std::dec << " Unknown Type" << std::endl;
00075 }
00076 }
00077 }
00078 }
00079
00080 std::vector<DetId> matrixECALIds(const DetId& det,int ieta,int iphi,
00081 const CaloGeometry* geo,
00082 const CaloTopology* caloTopology,
00083 bool debug) {
00084
00085 std::vector<DetId> vdets;
00086 spr::matrixECALIds(det, ieta, iphi, geo, caloTopology, vdets, debug);
00087 return vdets;
00088 }
00089
00090 std::vector<DetId> matrixECALIds(const DetId& det, double dR,
00091 const GlobalVector& trackMom,
00092 const CaloGeometry* geo,
00093 const CaloTopology* caloTopology,
00094 bool debug) {
00095
00096 GlobalPoint core;
00097 if (det.subdetId() == EcalEndcap) {
00098 EEDetId EEid = EEDetId(det);
00099 core = geo->getPosition(EEid);
00100 } else {
00101 EBDetId EBid = EBDetId(det);
00102 core = geo->getPosition(EBid);
00103 }
00104 int ietaphi = (int)(dR/2.0)+1;
00105 std::vector<DetId> vdets, vdetx;
00106 spr::matrixECALIds(det, ietaphi, ietaphi, geo, caloTopology, vdets, debug);
00107 for (unsigned int i=0; i<vdets.size(); ++i) {
00108 GlobalPoint rpoint;
00109 if (vdets[i].subdetId() == EcalEndcap) {
00110 EEDetId EEid = EEDetId(vdets[i]);
00111 rpoint = geo->getPosition(EEid);
00112 } else {
00113 EBDetId EBid = EBDetId(vdets[i]);
00114 rpoint = geo->getPosition(EBid);
00115 }
00116 if (spr::getDistInPlaneTrackDir(core, trackMom, rpoint)<dR) {
00117 vdetx.push_back(vdets[i]);
00118 }
00119 }
00120
00121 if (debug) {
00122 std::cout << "matrixECALIds::Final List of cells for dR " << dR
00123 << " is with " << vdetx.size() << " from original list of "
00124 << vdets.size() << std::endl;
00125 for (unsigned int i=0; i < vdetx.size(); ++i) {
00126 if (vdetx[i].subdetId() == EcalBarrel) {
00127 EBDetId id = vdetx[i];
00128 std::cout << "matrixECALIds::Cell " << i << " 0x" << std::hex
00129 << vdetx[i]() << std::dec << " " << id << std::endl;
00130 } else if (vdetx[i].subdetId() == EcalEndcap) {
00131 EEDetId id = vdetx[i];
00132 std::cout << "matrixECALIds::Cell " << i << " 0x" << std::hex
00133 << vdetx[i]() << std::dec<< " " << id << std::endl;
00134 } else {
00135 std::cout << "matrixECALIds::Cell " << i << " 0x" << std::hex
00136 << vdetx[i]() << std::dec << " Unknown Type" << std::endl;
00137 }
00138 }
00139 }
00140 return vdetx;
00141 }
00142
00143 void matrixECALIds(const DetId& det, int ietaE,int ietaW,int iphiN,int iphiS,
00144 const CaloGeometry* geo, const CaloTopology* caloTopology,
00145 std::vector<DetId>& vdets, bool debug) {
00146
00147 const CaloSubdetectorTopology *barrelTopo = (caloTopology->getSubdetectorTopology(DetId::Ecal,EcalBarrel));
00148 const CaloSubdetectorTopology *endcapTopo = (caloTopology->getSubdetectorTopology(DetId::Ecal,EcalEndcap));
00149 const EcalBarrelGeometry *barrelGeom = (dynamic_cast< const EcalBarrelGeometry *> (geo->getSubdetectorGeometry(DetId::Ecal,EcalBarrel)));
00150 const EcalEndcapGeometry *endcapGeom = (dynamic_cast< const EcalEndcapGeometry *> (geo->getSubdetectorGeometry(DetId::Ecal,EcalEndcap)));
00151
00152 if (debug) {
00153 std::cout << "matrixECALIds::Add " << ietaE << "|" << ietaW
00154 << " rows and " << iphiN << "|" << iphiS
00155 << " columns of cells for 1 cell" << std::endl;
00156 if (det.subdetId() == EcalBarrel) {
00157 EBDetId id = det;
00158 std::cout << "matrixECALIds::Cell 0x" << std::hex << det() << std::dec
00159 << " " << id << std::endl;
00160 } else if (det.subdetId() == EcalEndcap) {
00161 EEDetId id = det;
00162 std::cout << "matrixECALIds::Cell 0x" << std::hex << det() << std::dec
00163 << " " << id << std::endl;
00164 } else {
00165 std::cout << "matrixECALIds::Cell 0x" << std::hex << det() << std::dec
00166 << " Unknown Type" << std::endl;
00167 }
00168 }
00169
00170 std::vector<DetId> dets(1,det);
00171 std::vector<CaloDirection> dirs(1,NORTH);
00172 std::vector<int> jetaE(1,ietaE), jetaW(1,ietaW);
00173 std::vector<int> jphiN(1,iphiN), jphiS(1,iphiS);
00174 vdets = spr::newECALIdNS(dets, 0, jetaE, jetaW, jphiN, jphiS, dirs,
00175 *barrelTopo, *endcapTopo, *barrelGeom,
00176 *endcapGeom, debug);
00177 dirs[0] = SOUTH;
00178 std::vector<DetId> vdetS = spr::newECALIdNS(dets, 0, jetaE, jetaW, jphiN,
00179 jphiS, dirs, *barrelTopo,
00180 *endcapTopo, *barrelGeom,
00181 *endcapGeom, debug);
00182 for (unsigned int i1=0; i1<vdetS.size(); i1++) {
00183 if (std::count(vdets.begin(),vdets.end(),vdetS[i1]) == 0)
00184 vdets.push_back(vdetS[i1]);
00185 }
00186
00187 unsigned int ndet = (ietaE+ietaW+1)*(iphiN+iphiS+1);
00188 if (vdets.size() != ndet) {
00189 std::vector<DetId> vdetExtra;
00190 spr::extraIds(det, vdets, ietaE, ietaW, iphiN, iphiS,
00191 *barrelGeom, *endcapGeom, vdetExtra, debug);
00192 if (vdetExtra.size() > 0)
00193 vdets.insert(vdets.end(), vdetExtra.begin(), vdetExtra.end());
00194 }
00195
00196 if (debug) {
00197 std::cout << "matrixECALIds::Total number of cells found is "
00198 << vdets.size() << std::endl;
00199 for (unsigned int i1=0; i1<vdets.size(); i1++) {
00200 if (vdets[i1].subdetId() == EcalBarrel) {
00201 EBDetId id = vdets[i1];
00202 std::cout << "matrixECALIds::Cell " << i1 << " 0x" << std::hex
00203 << vdets[i1]() << std::dec << " " << id << std::endl;
00204 } else if (vdets[i1].subdetId() == EcalEndcap) {
00205 EEDetId id = vdets[i1];
00206 std::cout << "matrixECALIds::Cell " << i1 << " 0x" << std::hex
00207 << vdets[i1]() << std::dec << " " << id << std::endl;
00208 } else {
00209 std::cout << "matrixECALIds::Cell " << i1 << " 0x" << std::hex
00210 << vdets[i1]() << std::dec << " Unknown Type" << std::endl;
00211 }
00212 }
00213 }
00214 }
00215
00216 std::vector<DetId> matrixECALIds(const DetId& det, int ietaE, int ietaW,
00217 int iphiN, int iphiS,
00218 const CaloGeometry* geo,
00219 const CaloTopology* caloTopology,
00220 bool debug) {
00221
00222 std::vector<DetId> vdets;
00223 spr::matrixECALIds(det, ietaE, ietaW, iphiN, iphiS, geo, caloTopology,
00224 vdets, debug);
00225 return vdets;
00226 }
00227
00228 std::vector<DetId> newECALIdNS(std::vector<DetId>& dets, unsigned int last,
00229 int ieta, int iphi,
00230 std::vector<CaloDirection>& dir,
00231 const CaloSubdetectorTopology& barrelTopo,
00232 const CaloSubdetectorTopology& endcapTopo,
00233 const EcalBarrelGeometry& barrelGeom,
00234 const EcalEndcapGeometry& endcapGeom,
00235 bool debug) {
00236
00237 if (debug) {
00238 std::cout << "newECALIdNS::Add " << iphi << " columns of cells for "
00239 << (dets.size()-last) << " cells (last " << last << ")"
00240 << std::endl;
00241 for (unsigned int i1=last; i1<dets.size(); i1++) {
00242 if (dets[i1].subdetId() == EcalBarrel) {
00243 EBDetId id = dets[i1];
00244 std::cout << "newECALIdNS::Cell " << i1 << " " << id << " along "
00245 << dir[i1] << std::endl;
00246 } else if (dets[i1].subdetId() == EcalEndcap) {
00247 EEDetId id = dets[i1];
00248 std::cout << "newECALIdNS::Cell " << i1 << " " << id << " along "
00249 << dir[i1] << std::endl;
00250 } else {
00251 std::cout << "newECALIdNS::Cell " << i1 << " 0x" << std::hex
00252 << dets[i1]() << std::dec << " Unknown Type along "
00253 << dir[i1] << std::endl;
00254 }
00255 }
00256 }
00257
00258 std::vector<DetId> vdets;
00259 std::vector<CaloDirection> dirs;
00260 vdets.insert(vdets.end(), dets.begin(), dets.end());
00261 dirs.insert(dirs.end(), dir.begin(), dir.end());
00262
00263 std::vector<DetId> vdetE, vdetW;
00264 if (last == 0) {
00265 unsigned int ndet = vdets.size();
00266 std::vector<CaloDirection> dirE(ndet,EAST), dirW(ndet,WEST);
00267 vdetE = spr::newECALIdEW(dets, last, ieta, dirE, barrelTopo, endcapTopo,
00268 barrelGeom, endcapGeom, debug);
00269 vdetW = spr::newECALIdEW(dets, last, ieta, dirW, barrelTopo, endcapTopo,
00270 barrelGeom, endcapGeom, debug);
00271 for (unsigned int i1=0; i1<vdetW.size(); i1++) {
00272 if (std::count(vdets.begin(),vdets.end(),vdetW[i1]) == 0) {
00273 vdets.push_back(vdetW[i1]);
00274 dirs.push_back(dir[0]);
00275 }
00276 }
00277 for (unsigned int i1=0; i1<vdetE.size(); i1++) {
00278 if (std::count(vdets.begin(),vdets.end(),vdetE[i1]) == 0) {
00279 vdets.push_back(vdetE[i1]);
00280 dirs.push_back(dir[0]);
00281 }
00282 }
00283 if (debug) {
00284 std::cout <<"newECALIdNS::With Added cells along E/W results a set of "
00285 << (vdets.size()-dets.size()) << " new cells" << std::endl;
00286 for (unsigned int i1=dets.size(); i1<vdets.size(); i1++) {
00287 if (vdets[i1].subdetId() == EcalBarrel) {
00288 EBDetId id = vdets[i1];
00289 std::cout << "newECALIdNS::Cell " << i1 << " " << id << " along "
00290 << dirs[i1] << std::endl;
00291 } else if (vdets[i1].subdetId() == EcalEndcap) {
00292 EEDetId id = vdets[i1];
00293 std::cout << "newECALIdNS::Cell " << i1 << " " << id << " along "
00294 << dirs[i1] << std::endl;
00295 } else {
00296 std::cout << "newECALIdNS::Cell " << i1 << " 0x" << std::hex
00297 << vdets[i1]() << std::dec << " Unknown Type along "
00298 << dirs[i1] << std::endl;
00299 }
00300 }
00301 }
00302 }
00303
00304 unsigned int last0 = vdets.size();
00305 std::vector<DetId> vdetnew;
00306 std::vector<CaloDirection> dirnew;
00307 if (iphi > 0) {
00308 std::vector<DetId> vdetn(1);
00309 std::vector<CaloDirection> dirn(1);
00310 std::vector<CaloDirection> dirnE(1,EAST), dirnW(1,WEST);
00311 int flag=0;
00312 for (unsigned int i1=last; i1<dets.size(); i1++) {
00313 std::vector<DetId> cells;
00314 spr::simpleMove(dets[i1], dir[i1], barrelTopo, endcapTopo,
00315 barrelGeom, endcapGeom, cells, flag, debug);
00316
00317
00318
00319
00320
00321
00322
00323
00324
00325
00326
00327
00328
00329
00330
00331
00332
00333
00334
00335
00336
00337
00338
00339
00340
00341
00342
00343
00344 if (flag != 0) {
00345 if (std::count(vdets.begin(),vdets.end(),cells[0]) == 0) {
00346 vdetn[0] = cells[0];
00347 vdetnew.push_back(vdetn[0]);
00348 dirn[0] = dir[i1];
00349 if (flag < 0) {
00350 if (dirn[0] == NORTH) dirn[0] = SOUTH;
00351 else dirn[0] = NORTH;
00352 }
00353 dirnew.push_back(dirn[0]);
00354 vdetE = spr::newECALIdEW(vdetn, 0, ieta, dirnE, barrelTopo,
00355 endcapTopo, barrelGeom, endcapGeom,debug);
00356 vdetW = spr::newECALIdEW(vdetn, 0, ieta, dirnW, barrelTopo,
00357 endcapTopo, barrelGeom, endcapGeom,debug);
00358 for (unsigned int i2=0; i2<vdetW.size(); i2++) {
00359 if (std::count(vdets.begin(),vdets.end(),vdetW[i2]) == 0 &&
00360 std::count(vdetnew.begin(),vdetnew.end(),vdetW[i2]) == 0) {
00361 vdets.push_back(vdetW[i2]);
00362 dirs.push_back(dirn[0]);
00363 }
00364 }
00365 for (unsigned int i2=0; i2<vdetE.size(); i2++) {
00366 if (std::count(vdets.begin(),vdets.end(),vdetE[i2]) == 0 &&
00367 std::count(vdetnew.begin(),vdetnew.end(),vdetE[i2]) == 0) {
00368 vdets.push_back(vdetE[i2]);
00369 dirs.push_back(dirn[0]);
00370 }
00371 }
00372 }
00373 }
00374 }
00375 iphi--;
00376 last = vdets.size();
00377 for (unsigned int i2=0; i2<vdetnew.size(); i2++) {
00378 if (std::count(vdets.begin(),vdets.end(),vdetnew[i2]) == 0) {
00379 vdets.push_back(vdetnew[i2]);
00380 dirs.push_back(dirnew[i2]);
00381 }
00382 }
00383 if (debug) {
00384 std::cout << "newECALIdNS::Addition results a set of "
00385 << (vdets.size()-last0) << " new cells (last " << last0
00386 << ", iphi " << iphi << ")" << std::endl;
00387 for (unsigned int i1=last0; i1<vdets.size(); i1++) {
00388 if (vdets[i1].subdetId() == EcalBarrel) {
00389 EBDetId id = vdets[i1];
00390 std::cout << "newECALIdNS::Cell " << i1 << " " << id << " along "
00391 << dirs[i1] << std::endl;
00392 } else if (vdets[i1].subdetId() == EcalEndcap) {
00393 EEDetId id = vdets[i1];
00394 std::cout << "newECALIdNS::Cell " << i1 << " " << id << " along "
00395 << dirs[i1] << std::endl;
00396 } else {
00397 std::cout << "newECALIdNS::Cell " << i1 << " 0x" << std::hex
00398 << vdets[i1]() << std::dec << " Unknown Type along "
00399 << dirs[i1] << std::endl;
00400 }
00401 }
00402 }
00403 last0 = last;
00404 }
00405
00406 if (iphi > 0) {
00407 last = last0;
00408 return spr::newECALIdNS(vdets,last,ieta,iphi,dirs,barrelTopo,endcapTopo,barrelGeom,endcapGeom,debug);
00409 } else {
00410 if (debug) {
00411 std::cout << "newECALIdNS::Final list consists of " << vdets.size()
00412 << " cells" << std::endl;
00413 for (unsigned int i1=0; i1<vdets.size(); i1++) {
00414 if (vdets[i1].subdetId() == EcalBarrel) {
00415 EBDetId id = vdets[i1];
00416 std::cout << "newECALIdNS::Cell " << i1 << " " << id << std::endl;
00417 } else if (vdets[i1].subdetId() == EcalEndcap) {
00418 EEDetId id = vdets[i1];
00419 std::cout << "newECALIdNS::Cell " << i1 << " " << id << std::endl;
00420 } else {
00421 std::cout << "newECALIdNS::Cell " << i1 << " 0x" << std::hex
00422 << vdets[i1]() << std::dec << " Unknown Type"
00423 << std::endl;
00424 }
00425 }
00426 }
00427 return vdets;
00428 }
00429 }
00430
00431 std::vector<DetId> newECALIdNS(std::vector<DetId>& dets, unsigned int last,
00432 std::vector<int>& ietaE,
00433 std::vector<int>& ietaW,
00434 std::vector<int>& iphiN,
00435 std::vector<int>& iphiS,
00436 std::vector<CaloDirection>& dir,
00437 const CaloSubdetectorTopology& barrelTopo,
00438 const CaloSubdetectorTopology& endcapTopo,
00439 const EcalBarrelGeometry& barrelGeom,
00440 const EcalEndcapGeometry& endcapGeom,
00441 bool debug) {
00442
00443 if (debug) {
00444 std::cout << "newECALIdNS::Add columns of cells for "
00445 << (dets.size()-last) << " cells (last) " << last << std::endl;
00446 for (unsigned int i1=last; i1<dets.size(); i1++) {
00447 if (dets[i1].subdetId() == EcalBarrel) {
00448 EBDetId id = dets[i1];
00449 std::cout << "newECALIdNS::Cell " << i1 << " " << id << " along "
00450 << dir[i1] << " # " << iphiN[i1] << "|" << iphiS[i1]
00451 << std::endl;
00452 } else if (dets[i1].subdetId() == EcalEndcap) {
00453 EEDetId id = dets[i1];
00454 std::cout << "newECALIdNS::Cell " << i1 << " " << id << " along "
00455 << dir[i1] << " # " << iphiN[i1] << "|" << iphiS[i1]
00456 << std::endl;
00457 } else {
00458 std::cout << "newECALIdNS::Cell " << i1 << " 0x" << std::hex
00459 << dets[i1]() << std::dec << " Unknown Type along "
00460 << dir[i1] << " # " << iphiN[i1] << "|" << iphiS[i1]
00461 << std::endl;
00462 }
00463 }
00464 }
00465
00466 std::vector<DetId> vdets;
00467 std::vector<CaloDirection> dirs;
00468 std::vector<int> jetaE, jetaW, jphiN, jphiS;
00469 vdets.insert(vdets.end(), dets.begin(), dets.end());
00470 dirs.insert(dirs.end(), dir.begin(), dir.end());
00471 jetaE.insert(jetaE.end(), ietaE.begin(), ietaE.end());
00472 jetaW.insert(jetaW.end(), ietaW.begin(), ietaW.end());
00473 jphiN.insert(jphiN.end(), iphiN.begin(), iphiN.end());
00474 jphiS.insert(jphiS.end(), iphiS.begin(), iphiS.end());
00475 std::vector<DetId> vdetE, vdetW;
00476 if (last == 0) {
00477 unsigned int ndet = vdets.size();
00478 std::vector<CaloDirection> dirE(ndet,EAST), dirW(ndet,WEST);
00479 vdetE = spr::newECALIdEW(dets, last, ietaE, ietaW, dirE, barrelTopo,
00480 endcapTopo, barrelGeom, endcapGeom, debug);
00481 vdetW = spr::newECALIdEW(dets, last, ietaE, ietaW, dirW, barrelTopo,
00482 endcapTopo, barrelGeom, endcapGeom, debug);
00483 for (unsigned int i1=0; i1<vdetW.size(); i1++) {
00484 if (std::count(vdets.begin(),vdets.end(),vdetW[i1]) == 0) {
00485 vdets.push_back(vdetW[i1]);
00486 dirs.push_back(dir[0]);
00487 jetaE.push_back(0);
00488 jetaW.push_back(0);
00489 jphiN.push_back(iphiN[0]);
00490 jphiS.push_back(iphiS[0]);
00491 }
00492 }
00493 for (unsigned int i1=0; i1<vdetE.size(); i1++) {
00494 if (std::count(vdets.begin(),vdets.end(),vdetE[i1]) == 0) {
00495 vdets.push_back(vdetE[i1]);
00496 dirs.push_back(dir[0]);
00497 jetaE.push_back(0);
00498 jetaW.push_back(0);
00499 jphiN.push_back(iphiN[0]);
00500 jphiS.push_back(iphiS[0]);
00501 }
00502 }
00503 if (debug) {
00504 std::cout <<"newECALIdNS::With Added cells along E/W results a set of "
00505 << (vdets.size()-dets.size()) << " new cells" << std::endl;
00506 for (unsigned int i1=dets.size(); i1<vdets.size(); i1++) {
00507 if (vdets[i1].subdetId() == EcalBarrel) {
00508 EBDetId id = vdets[i1];
00509 std::cout << "newECALIdNS::Cell " << i1 << " " << id << " along "
00510 << dirs[i1] << std::endl;
00511 } else if (vdets[i1].subdetId() == EcalEndcap) {
00512 EEDetId id = vdets[i1];
00513 std::cout << "newECALIdNS::Cell " << i1 << " " << id << " along "
00514 << dirs[i1] << std::endl;
00515 } else {
00516 std::cout << "newECALIdNS::Cell " << i1 << " 0x" << std::hex
00517 << vdets[i1]() << std::dec << " Unknown Type along "
00518 << dirs[i1] << std::endl;
00519 }
00520 }
00521 }
00522 }
00523
00524 unsigned int last0 = vdets.size();
00525 std::vector<DetId> vdetnew;
00526 std::vector<CaloDirection> dirnew;
00527 std::vector<int> kphiN, kphiS, ketaE, ketaW;
00528 int kphi = 0;
00529 for (unsigned int i1=last; i1<dets.size(); i1++) {
00530 int iphi = iphiS[i1];
00531 if (dir[i1] == NORTH) iphi = iphiN[i1];
00532 if (iphi > 0) {
00533 std::vector<DetId> vdetn(1);
00534 std::vector<CaloDirection> dirn(1);
00535 std::vector<CaloDirection> dirnE(1,EAST), dirnW(1,WEST);
00536 int flag=0;
00537 std::vector<DetId> cells;
00538 spr::simpleMove(dets[i1], dir[i1], barrelTopo, endcapTopo,
00539 barrelGeom, endcapGeom, cells, flag, debug);
00540 iphi--;
00541 if (iphi > kphi) kphi = iphi;
00542 if (dir[i1] == NORTH) jphiN[i1] = iphi;
00543 else jphiS[i1] = iphi;
00544 if (flag != 0) {
00545 if (std::count(vdets.begin(),vdets.end(),cells[0]) == 0) {
00546 int kfiN = iphiN[i1];
00547 int kfiS = iphiS[i1];
00548 vdetn[0] = cells[0];
00549 vdetnew.push_back(vdetn[0]);
00550 dirn[0] = dir[i1];
00551 if (dir[i1] == NORTH) kfiN = iphi;
00552 else kfiS = iphi;
00553 if (flag < 0) {
00554 int ktmp = kfiS; kfiS = kfiN; kfiN = ktmp;
00555 if (dirn[0] == NORTH) dirn[0] = SOUTH;
00556 else dirn[0] = NORTH;
00557 }
00558 dirnew.push_back(dirn[0]);
00559 kphiN.push_back(kfiN); ketaE.push_back(ietaE[i1]);
00560 kphiS.push_back(kfiS); ketaW.push_back(ietaW[i1]);
00561 std::vector<int> ietE(1,ietaE[i1]), ietW(1,ietaW[i1]);
00562 vdetE = spr::newECALIdEW(vdetn, 0, ietE, ietW, dirnE, barrelTopo,
00563 endcapTopo, barrelGeom, endcapGeom,debug);
00564 vdetW = spr::newECALIdEW(vdetn, 0, ietE, ietW, dirnW, barrelTopo,
00565 endcapTopo, barrelGeom, endcapGeom,debug);
00566 for (unsigned int i2=0; i2<vdetW.size(); i2++) {
00567 if (std::count(vdets.begin(),vdets.end(),vdetW[i2]) == 0 &&
00568 std::count(vdetnew.begin(),vdetnew.end(),vdetW[i2]) == 0) {
00569 vdets.push_back(vdetW[i2]);
00570 dirs.push_back(dirn[0]);
00571 jetaE.push_back(0); jphiN.push_back(kfiN);
00572 jetaW.push_back(0); jphiS.push_back(kfiS);
00573 }
00574 }
00575 for (unsigned int i2=0; i2<vdetE.size(); i2++) {
00576 if (std::count(vdets.begin(),vdets.end(),vdetE[i2]) == 0 &&
00577 std::count(vdetnew.begin(),vdetnew.end(),vdetE[i2]) == 0) {
00578 vdets.push_back(vdetE[i2]);
00579 dirs.push_back(dirn[0]);
00580 jetaE.push_back(0); jphiN.push_back(kfiN);
00581 jetaW.push_back(0); jphiS.push_back(kfiS);
00582 }
00583 }
00584 }
00585 }
00586 }
00587 }
00588 last = vdets.size();
00589 for (unsigned int i2=0; i2<vdetnew.size(); i2++) {
00590 if (std::count(vdets.begin(),vdets.end(),vdetnew[i2]) == 0) {
00591 vdets.push_back(vdetnew[i2]);
00592 dirs.push_back(dirnew[i2]);
00593 jetaE.push_back(ketaE[i2]);
00594 jetaW.push_back(ketaW[i2]);
00595 jphiN.push_back(kphiN[i2]);
00596 jphiS.push_back(kphiS[i2]);
00597 }
00598 }
00599 if (debug) {
00600 std::cout << "newECALIdNS::Addition results a set of "
00601 << (vdets.size()-last0) << " new cells (last " << last0
00602 << ", iphi " << kphi << ")" << std::endl;
00603 for (unsigned int i1=last0; i1<vdets.size(); i1++) {
00604 if (vdets[i1].subdetId() == EcalBarrel) {
00605 EBDetId id = vdets[i1];
00606 std::cout << "newECALIdNS::Cell " << i1 << " " << id << " along "
00607 << dirs[i1] << " iphi " << jphiN[i1] << "|" << jphiS[i1]
00608 << " ieta " << jetaE[i1] << "|" << jetaW[i1] << std::endl;
00609 } else if (vdets[i1].subdetId() == EcalEndcap) {
00610 EEDetId id = vdets[i1];
00611 std::cout << "newECALIdNS::Cell " << i1 << " " << id << " along "
00612 << dirs[i1] << " iphi " << jphiN[i1] << "|" << jphiS[i1]
00613 << " ieta " << jetaE[i1] << "|" << jetaW[i1] << std::endl;
00614 } else {
00615 std::cout << "newECALIdNS::Cell " << i1 << " 0x" << std::hex
00616 << vdets[i1]() << std::dec << " Unknown Type along "
00617 << dirs[i1] << " iphi " << jphiN[i1] << "|" << jphiS[i1]
00618 << " ieta " << jetaE[i1] << "|" << jetaW[i1] << std::endl;
00619 }
00620 }
00621 }
00622 last0 = last;
00623
00624 if (kphi > 0) {
00625 last = last0;
00626 return spr::newECALIdNS(vdets, last, jetaE, jetaW, jphiN, jphiS, dirs,
00627 barrelTopo, endcapTopo, barrelGeom, endcapGeom,
00628 debug);
00629 } else {
00630 if (debug) {
00631 std::cout << "newECALIdNS::Final list consists of " << vdets.size()
00632 << " cells" << std::endl;
00633 for (unsigned int i1=0; i1<vdets.size(); i1++) {
00634 if (vdets[i1].subdetId() == EcalBarrel) {
00635 EBDetId id = vdets[i1];
00636 std::cout << "newECALIdNS::Cell " << i1 << " " << id << std::endl;
00637 } else if (vdets[i1].subdetId() == EcalEndcap) {
00638 EEDetId id = vdets[i1];
00639 std::cout << "newECALIdNS::Cell " << i1 << " " << id << std::endl;
00640 } else {
00641 std::cout << "newECALIdNS::Cell " << i1 << " 0x" << std::hex
00642 << vdets[i1]() << std::dec << " Unknown Type"
00643 << std::endl;
00644 }
00645 }
00646 }
00647 return vdets;
00648 }
00649 }
00650
00651 std::vector<DetId> newECALIdEW(std::vector<DetId>& dets, unsigned int last,
00652 int ieta, std::vector<CaloDirection>& dir,
00653 const CaloSubdetectorTopology& barrelTopo,
00654 const CaloSubdetectorTopology& endcapTopo,
00655 const EcalBarrelGeometry& barrelGeom,
00656 const EcalEndcapGeometry& endcapGeom,
00657 bool debug) {
00658
00659 if (debug) {
00660 std::cout << "newECALIdEW::Add " << ieta << " rows of cells for "
00661 << last << ":" << dets.size() << ":" << (dets.size()-last)
00662 << " cells" << std::endl;
00663 for (unsigned int i1=last; i1<dets.size(); i1++) {
00664 if (dets[i1].subdetId() == EcalBarrel) {
00665 EBDetId id = dets[i1];
00666 std::cout << "newECALIdEW::Cell " << i1 << " " << id << " along "
00667 << dir[i1] << std::endl;
00668 } else if (dets[i1].subdetId() == EcalEndcap) {
00669 EEDetId id = dets[i1];
00670 std::cout << "newECALIdEW::Cell " << i1 << " " << id << " along "
00671 << dir[i1] << std::endl;
00672 } else {
00673 std::cout << "newECALIdEW::Cell " << i1 << " 0x" << std::hex
00674 << dets[i1]() << std::dec << " Unknown Type along "
00675 << dir[i1] << std::endl;
00676 }
00677 }
00678 }
00679
00680 std::vector<DetId> vdets; vdets.clear();
00681 std::vector<CaloDirection> dirs; dirs.clear();
00682 vdets.insert(vdets.end(), dets.begin(), dets.end());
00683 dirs.insert(dirs.end(), dir.begin(), dir.end());
00684
00685 if (ieta > 0) {
00686 for (unsigned int i1=last; i1<dets.size(); i1++) {
00687 int flag = 0;
00688 std::vector<DetId> cells;
00689 spr::simpleMove(dets[i1], dir[i1], barrelTopo, endcapTopo,
00690 barrelGeom, endcapGeom, cells, flag, debug);
00691
00692
00693
00694
00695
00696
00697
00698
00699
00700
00701
00702
00703
00704
00705
00706
00707
00708
00709
00710
00711
00712
00713
00714
00715
00716
00717
00718
00719 if (flag != 0) {
00720 if (std::count(vdets.begin(),vdets.end(),cells[0]) == 0) {
00721 CaloDirection dirn = dir[i1];
00722 if (flag < 0) {
00723 if (dirn == EAST) dirn = WEST;
00724 else dirn = EAST;
00725 }
00726 vdets.push_back(cells[0]);
00727 dirs.push_back(dirn);
00728 }
00729 }
00730 }
00731 ieta--;
00732 }
00733
00734 if (debug) {
00735 std::cout << "newECALIdEW::Addition results a set of "
00736 << (vdets.size()-dets.size()) << " new cells" << std::endl;
00737 for (unsigned int i1=dets.size(); i1<vdets.size(); i1++) {
00738 if (vdets[i1].subdetId() == EcalBarrel) {
00739 EBDetId id = vdets[i1];
00740 std::cout << "newECALIdEW::Cell " << i1 << " " << id << " along "
00741 << dirs[i1] << std::endl;
00742 } else if (vdets[i1].subdetId() == EcalEndcap) {
00743 EEDetId id = vdets[i1];
00744 std::cout << "newECALIdEW::Cell " << i1 << " " << id << " along "
00745 << dirs[i1] << std::endl;
00746 } else {
00747 std::cout << "newECALIdEW::Cell " << i1 << " 0x" << std::hex
00748 << vdets[i1]() << std::dec << " Unknown Type along "
00749 << dirs[i1] << std::endl;
00750 }
00751 }
00752 }
00753
00754 if (ieta > 0) {
00755 last = dets.size();
00756 return spr::newECALIdEW(vdets,last,ieta,dirs,barrelTopo,endcapTopo,barrelGeom,endcapGeom,debug);
00757 } else {
00758 if (debug) {
00759 std::cout << "newECALIdEW::Final list (EW) consists of " <<vdets.size()
00760 << " cells" << std::endl;
00761 for (unsigned int i1=0; i1<vdets.size(); i1++) {
00762 if (vdets[i1].subdetId() == EcalBarrel) {
00763 EBDetId id = vdets[i1];
00764 std::cout << "newECALIdEW::Cell " << i1 << " " << id << std::endl;
00765 } else if (vdets[i1].subdetId() == EcalEndcap) {
00766 EEDetId id = vdets[i1];
00767 std::cout << "newECALIdEW::Cell " << i1 << " " << id << std::endl;
00768 } else {
00769 std::cout << "newECALIdEW::Cell " << i1 << " 0x" << std::hex
00770 << vdets[i1]() <<std::dec << " Unknown Type" <<std::endl;
00771 }
00772 }
00773 }
00774 return vdets;
00775 }
00776 }
00777
00778 std::vector<DetId> newECALIdEW(std::vector<DetId>& dets, unsigned int last,
00779 std::vector<int>& ietaE,
00780 std::vector<int>& ietaW,
00781 std::vector<CaloDirection>& dir,
00782 const CaloSubdetectorTopology& barrelTopo,
00783 const CaloSubdetectorTopology& endcapTopo,
00784 const EcalBarrelGeometry& barrelGeom,
00785 const EcalEndcapGeometry& endcapGeom,
00786 bool debug) {
00787
00788 if (debug) {
00789 std::cout << "newECALIdEW::Add " << ietaE[0] << "|" << ietaW[0]
00790 << " rows of cells for " << (dets.size()-last)
00791 << " cells (last " << last << ")" << std::endl;
00792 for (unsigned int i1=last; i1<dets.size(); i1++) {
00793 if (dets[i1].subdetId() == EcalBarrel) {
00794 EBDetId id = dets[i1];
00795 std::cout << "newECALIdEW::Cell " << i1 << " " << id << " along "
00796 << dir[i1] << std::endl;
00797 } else if (dets[i1].subdetId() == EcalEndcap) {
00798 EEDetId id = dets[i1];
00799 std::cout << "newECALIdEW::Cell " << i1 << " " << id << " along "
00800 << dir[i1] << std::endl;
00801 } else {
00802 std::cout << "newECALIdEW::Cell " << i1 << " 0x" << std::hex
00803 << dets[i1]() << std::dec << " Unknown Type along "
00804 << dir[i1] << std::endl;
00805 }
00806 }
00807 }
00808
00809 std::vector<DetId> vdets;
00810 vdets.insert(vdets.end(), dets.begin(), dets.end());
00811 std::vector<CaloDirection> dirs;
00812 dirs.insert(dirs.end(), dir.begin(), dir.end());
00813 std::vector<int> jetaE, jetaW;
00814 jetaE.insert(jetaE.end(), ietaE.begin(), ietaE.end());
00815 jetaW.insert(jetaW.end(), ietaW.begin(), ietaW.end());
00816 int keta = 0;
00817 for (unsigned int i1=last; i1<dets.size(); i1++) {
00818 int ieta = ietaW[i1];
00819 if (dir[i1] == EAST) ieta = ietaE[i1];
00820 if (ieta > 0) {
00821 int flag=0;
00822 std::vector<DetId> cells;
00823 spr::simpleMove(dets[i1], dir[i1], barrelTopo, endcapTopo,
00824 barrelGeom, endcapGeom, cells, flag, debug);
00825 ieta--;
00826 if (ieta > keta) keta = ieta;
00827 if (dir[i1] == EAST) jetaE[i1] = ieta;
00828 else jetaW[i1] = ieta;
00829 if (flag != 0) {
00830 if (std::count(vdets.begin(),vdets.end(),cells[0]) == 0) {
00831 vdets.push_back(cells[0]);
00832 CaloDirection dirn = dir[i1];
00833 int ketaE = ietaE[i1];
00834 int ketaW = ietaW[i1];
00835 if (dirn == EAST) ketaE = ieta;
00836 else ketaW = ieta;
00837 if (flag < 0) {
00838 int ktmp = ketaW; ketaW = ketaE; ketaE = ktmp;
00839 if (dirn == EAST) dirn = WEST;
00840 else dirn = EAST;
00841 }
00842 dirs.push_back(dirn);
00843 jetaE.push_back(ketaE);
00844 jetaW.push_back(ketaW);
00845 }
00846 }
00847 }
00848 }
00849
00850 if (debug) {
00851 std::cout << "newECALIdEW::Addition results a set of "
00852 << (vdets.size()-dets.size()) << " new cells (last "
00853 << dets.size() << ", ieta " << keta << ")" << std::endl;
00854 for (unsigned int i1=dets.size(); i1<vdets.size(); i1++) {
00855 if (vdets[i1].subdetId() == EcalBarrel) {
00856 EBDetId id = vdets[i1];
00857 std::cout << "newECALIdEW::Cell " << i1 << " " << id << std::endl;
00858 } else if (vdets[i1].subdetId() == EcalEndcap) {
00859 EEDetId id = vdets[i1];
00860 std::cout << "newECALIdEW::Cell " << i1 << " " << id << std::endl;
00861 } else {
00862 std::cout << "newECALIdEW::Cell " << i1 << " 0x" << std::hex
00863 << vdets[i1]() << std::dec << " Unknown Type" << std::endl;
00864 }
00865 }
00866 }
00867
00868 if (keta > 0) {
00869 last = dets.size();
00870 return spr::newECALIdEW(vdets, last, jetaE, jetaW, dirs, barrelTopo,
00871 endcapTopo, barrelGeom, endcapGeom, debug);
00872 } else {
00873 if (debug) {
00874 std::cout << "newECALIdEW::Final list (EW) consists of " <<vdets.size()
00875 << " cells" << std::endl;
00876 for (unsigned int i1=0; i1<vdets.size(); i1++) {
00877 if (vdets[i1].subdetId() == EcalBarrel) {
00878 EBDetId id = vdets[i1];
00879 std::cout << "newECALIdEW::Cell " << i1 << " " << id << std::endl;
00880 } else if (vdets[i1].subdetId() == EcalEndcap) {
00881 EEDetId id = vdets[i1];
00882 std::cout << "newECALIdEW::Cell " << i1 << " " << id << std::endl;
00883 } else {
00884 std::cout << "newECALIdEW::Cell " << i1 << " 0x" << std::hex
00885 << vdets[i1]() <<std::dec << " Unknown Type" <<std::endl;
00886 }
00887 }
00888 }
00889 return vdets;
00890 }
00891 }
00892
00893 void simpleMove(DetId& det, const CaloDirection& dir,
00894 const CaloSubdetectorTopology& barrelTopo,
00895 const CaloSubdetectorTopology& endcapTopo,
00896 const EcalBarrelGeometry& barrelGeom,
00897 const EcalEndcapGeometry& endcapGeom,
00898 std::vector<DetId>& cells, int& ok, bool debug) {
00899
00900 DetId cell;
00901 ok = 0;
00902 if (det.subdetId() == EcalBarrel) {
00903 EBDetId detId = det;
00904 std::vector<DetId> neighbours = barrelTopo.getNeighbours(detId,dir);
00905 if (neighbours.size()>0 && !neighbours[0].null()) {
00906 cells.push_back(neighbours[0]);
00907 cell = neighbours[0];
00908 ok = 1;
00909 } else {
00910 const int ietaAbs ( detId.ietaAbs() ) ;
00911 if (EBDetId::MAX_IETA == ietaAbs ) {
00912
00913 const EcalBarrelGeometry::OrderedListOfEEDetId&
00914 ol( * barrelGeom.getClosestEndcapCells(detId) ) ;
00915
00916 cell = *(ol.begin() );
00917 neighbours = endcapTopo.getNeighbours(cell,dir);
00918 if (neighbours.size()>0 && !neighbours[0].null()) ok = 1;
00919 else ok =-1;
00920 for (EcalBarrelGeometry::OrderedListOfEEDetId::const_iterator iptr=ol.begin(); iptr != ol.end(); ++iptr)
00921 cells.push_back(*iptr);
00922 }
00923 }
00924 } else if (det.subdetId() == EcalEndcap) {
00925 EEDetId detId = det;
00926 std::vector<DetId> neighbours = endcapTopo.getNeighbours(detId,dir);
00927 if (neighbours.size()>0 && !neighbours[0].null()) {
00928 cells.push_back(neighbours[0]);
00929 cell = neighbours[0];
00930 ok = 1;
00931 } else {
00932
00933 const int iphi ( detId.iPhiOuterRing() ) ;
00934 if (iphi!= 0) {
00935
00936 const EcalEndcapGeometry::OrderedListOfEBDetId&
00937 ol( * endcapGeom.getClosestBarrelCells(detId) ) ;
00938
00939 cell = *(ol.begin() );
00940 neighbours = barrelTopo.getNeighbours(cell,dir);
00941 if (neighbours.size()>0 && !neighbours[0].null()) ok = 1;
00942 else ok =-1;
00943 for (EcalEndcapGeometry::OrderedListOfEBDetId::const_iterator iptr=ol.begin(); iptr != ol.end(); ++iptr)
00944 cells.push_back(*iptr);
00945 }
00946 }
00947 }
00948 if (debug) {
00949 std::cout << "simpleMove:: Move DetId 0x" << std::hex << det()
00950 << std::dec << " along " << dir << " to get 0x" << std::hex
00951 << cell() << std::dec << " with flag " << ok << " # "
00952 << cells.size();
00953 for (unsigned int i1=0; i1<cells.size(); ++i1)
00954 std::cout << " " << std::hex << cells[0]() << std::dec;
00955 std::cout << std::endl;
00956 }
00957 }
00958
00959 void extraIds(const DetId& det, std::vector<DetId>& dets, int ietaE, int ietaW, int iphiN, int iphiS, const EcalBarrelGeometry& barrelGeom, const EcalEndcapGeometry& endcapGeom, std::vector<DetId>& cells, bool debug) {
00960
00961 if (det.subdetId() == EcalBarrel) {
00962 EBDetId id = det;
00963 if (debug) std::cout << "extraIds::Cell " << id << " rows " << ietaW
00964 << "|" << ietaE << " columns " << iphiS << "|"
00965 << iphiN << std::endl;
00966 int etaC = id.ietaAbs();
00967 int phiC = id.iphi();
00968 int zsid = id.zside();
00969 for (int eta = -ietaW; eta <= ietaE; ++eta) {
00970 for (int phi = -iphiS; phi <= iphiN; ++phi) {
00971 int iphi = phiC+phi;
00972 if (iphi < 0) iphi += 360;
00973 else if (iphi > 360) iphi -= 360;
00974 int ieta = zsid*(etaC+eta);
00975 if (EBDetId::validDetId(ieta,iphi)) {
00976 id = EBDetId(ieta,iphi);
00977 if (barrelGeom.present(id)) {
00978 if (std::count(dets.begin(),dets.end(),(DetId)id) == 0) {
00979 cells.push_back((DetId)id);
00980 }
00981 }
00982 }
00983 }
00984 }
00985 } else if (det.subdetId() == EcalEndcap) {
00986 EEDetId id = det;
00987 if (debug) std::cout << "extraIds::Cell " << id << " rows " << ietaW
00988 << "|" << ietaE << " columns " << iphiS << "|"
00989 << iphiN << std::endl;
00990 int ixC = id.ix();
00991 int iyC = id.iy();
00992 int zsid = id.zside();
00993 for (int kx = -ietaW; kx <= ietaE; ++kx) {
00994 for (int ky = -iphiS; ky <= iphiN; ++ky) {
00995 int ix = ixC+kx;
00996 int iy = iyC+ky;
00997 if (EEDetId::validDetId(ix,iy,zsid)) {
00998 id = EEDetId(ix,iy,zsid);
00999 if (endcapGeom.present(id)) {
01000 if (std::count(dets.begin(),dets.end(),(DetId)id) == 0) {
01001 cells.push_back((DetId)id);
01002 }
01003 }
01004 }
01005 }
01006 }
01007 }
01008
01009 if (debug) {
01010 std::cout << "extraIds:: finds " << cells.size() << " new cells"
01011 << std::endl;
01012 for (unsigned int i1=0; i1<cells.size(); ++i1) {
01013 if (cells[i1].subdetId() == EcalBarrel) {
01014 EBDetId id = cells[i1];
01015 std::cout << "extraIds::Cell " << i1 << " " << id << std::endl;
01016 } else if (cells[i1].subdetId() == EcalEndcap) {
01017 EEDetId id = cells[i1];
01018 std::cout << "ectraIds::Cell " << i1 << " " << id << std::endl;
01019 } else {
01020 std::cout << "extraIds::Cell " << i1 << " 0x" << std::hex
01021 << cells[i1]() <<std::dec << " Unknown Type" <<std::endl;
01022 }
01023 }
01024 }
01025 }
01026 }
1.7.3