File indexing completed on 2024-04-06 11:59:10
0001 #include "Calibration/IsolatedParticles/interface/DebugInfo.h"
0002 #include "Calibration/IsolatedParticles/interface/FindCaloHit.h"
0003 #include "Calibration/IsolatedParticles/interface/FindEtaPhi.h"
0004 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0005 #include <sstream>
0006
0007 namespace spr {
0008
0009 template <typename T>
0010 double eECALmatrix(const DetId& detId,
0011 edm::Handle<T>& hitsEB,
0012 edm::Handle<T>& hitsEE,
0013 const CaloGeometry* geo,
0014 const CaloTopology* caloTopology,
0015 int ieta,
0016 int iphi,
0017 double ebThr,
0018 double eeThr,
0019 double tMin,
0020 double tMax,
0021 bool debug) {
0022 std::vector<DetId> vdets;
0023 spr::matrixECALIds(detId, ieta, iphi, geo, caloTopology, vdets, false);
0024
0025 if (debug) {
0026 edm::LogVerbatim("IsoTrack") << "Inside eECALmatrix " << 2 * ieta + 1 << "X" << 2 * iphi + 1 << " nXtals "
0027 << vdets.size();
0028 spr::debugEcalDets(0, vdets);
0029 }
0030
0031 return spr::energyECAL(vdets, hitsEB, hitsEE, ebThr, eeThr, tMin, tMax, debug);
0032 }
0033
0034 template <typename T>
0035 double eECALmatrix(const DetId& detId,
0036 edm::Handle<T>& hitsEB,
0037 edm::Handle<T>& hitsEE,
0038 const CaloGeometry* geo,
0039 const CaloTopology* caloTopology,
0040 const EcalTrigTowerConstituentsMap& ttMap,
0041 int ieta,
0042 int iphi,
0043 double ebThr,
0044 double eeThr,
0045 double tMin,
0046 double tMax,
0047 bool debug) {
0048 std::vector<DetId> vdets;
0049 spr::matrixECALIds(detId, ieta, iphi, geo, caloTopology, vdets, debug);
0050
0051 if (debug) {
0052 edm::LogVerbatim("IsoTrack") << "Inside eECALmatrix " << 2 * ieta + 1 << "X" << 2 * iphi + 1 << " nXtals "
0053 << vdets.size();
0054 }
0055
0056 return spr::energyECAL(vdets, hitsEB, hitsEE, ttMap, ebThr, eeThr, tMin, tMax, debug);
0057 }
0058
0059 template <typename T>
0060 double eECALmatrix(const DetId& detId,
0061 edm::Handle<T>& hitsEB,
0062 edm::Handle<T>& hitsEE,
0063 const CaloGeometry* geo,
0064 const CaloTopology* caloTopology,
0065 int ietaE,
0066 int ietaW,
0067 int iphiN,
0068 int iphiS,
0069 double ebThr,
0070 double eeThr,
0071 double tMin,
0072 double tMax,
0073 bool debug) {
0074 std::vector<DetId> vdets;
0075 spr::matrixECALIds(detId, ietaE, ietaW, iphiN, iphiS, geo, caloTopology, vdets, debug);
0076
0077 if (debug) {
0078 edm::LogVerbatim("IsoTrack") << "Inside eECALmatrix " << ietaE + ietaW + 1 << "X" << iphiN + iphiS + 1
0079 << " nXtals " << vdets.size();
0080 }
0081
0082 return spr::energyECAL(vdets, hitsEB, hitsEE, ebThr, eeThr, tMin, tMax, debug);
0083 }
0084
0085 template <typename T>
0086 void hitECALmatrix(CaloNavigator<DetId>& navigator,
0087 edm::Handle<T>& hits,
0088 int ieta,
0089 int iphi,
0090 std::vector<typename T::const_iterator>& hitlist,
0091 bool debug) {
0092 DetId thisDet;
0093
0094 for (int dx = -ieta; dx < ieta + 1; ++dx) {
0095 for (int dy = -iphi; dy < iphi + 1; ++dy) {
0096
0097 thisDet = navigator.offsetBy(dx, dy);
0098
0099
0100 navigator.home();
0101
0102 if (thisDet != DetId(0)) {
0103 std::vector<typename T::const_iterator> hit;
0104 spr::findHit(hits, thisDet, hit, debug);
0105 std::ostringstream st1;
0106 if (debug && !hit.empty()) {
0107 if (thisDet.subdetId() == EcalBarrel) {
0108 EBDetId id = thisDet;
0109 st1 << "hitECALmatrix::Cell 0x" << std::hex << thisDet() << std::dec << " " << id;
0110 } else if (thisDet.subdetId() == EcalEndcap) {
0111 EEDetId id = thisDet;
0112 st1 << "hitECALmatrix::Cell 0x" << std::hex << thisDet() << std::dec << " " << id;
0113 } else {
0114 st1 << "hitECALMatrix::Cell 0x" << std::hex << thisDet() << std::dec << " Unknown Type";
0115 }
0116 }
0117 for (unsigned int ihit = 0; ihit < hit.size(); ihit++) {
0118 if (hit[ihit] != hits->end()) {
0119 hitlist.push_back(hit[ihit]);
0120 if (debug)
0121 st1 << " hit " << ihit << " " << hit[ihit]->energy();
0122 }
0123 }
0124 if (debug && !hit.empty())
0125 edm::LogVerbatim("IsoTrack") << st1.str();
0126 }
0127
0128 }
0129 }
0130 }
0131
0132 template <typename T>
0133 double energyECAL(std::vector<DetId>& vdets,
0134 edm::Handle<T>& hitsEB,
0135 edm::Handle<T>& hitsEE,
0136 double ebThr,
0137 double eeThr,
0138 double tMin,
0139 double tMax,
0140 bool debug) {
0141 double energySum = 0.0;
0142 for (unsigned int i1 = 0; i1 < vdets.size(); i1++) {
0143 if (vdets[i1] != DetId(0)) {
0144 std::vector<typename T::const_iterator> hit;
0145 if (vdets[i1].subdetId() == EcalBarrel) {
0146 spr::findHit(hitsEB, vdets[i1], hit, debug);
0147 } else if (vdets[i1].subdetId() == EcalEndcap) {
0148 spr::findHit(hitsEE, vdets[i1], hit, debug);
0149 }
0150 std::ostringstream st1;
0151 if (debug)
0152 st1 << "Crystal 0x" << std::hex << vdets[i1]() << std::dec;
0153 double ener = 0, ethr = ebThr;
0154 if (vdets[i1].subdetId() != EcalBarrel)
0155 ethr = eeThr;
0156 for (unsigned int ihit = 0; ihit < hit.size(); ihit++) {
0157 double en = 0, tt = 0;
0158 if (vdets[i1].subdetId() == EcalBarrel) {
0159 if (hit[ihit] != hitsEB->end()) {
0160 en = hit[ihit]->energy();
0161 tt = hit[ihit]->time();
0162 }
0163 } else if (vdets[i1].subdetId() == EcalEndcap) {
0164 if (hit[ihit] != hitsEE->end()) {
0165 en = hit[ihit]->energy();
0166 tt = hit[ihit]->time();
0167 }
0168 }
0169 if (debug)
0170 st1 << " " << ihit << " " << en << " Thr " << ethr;
0171 if (tt > tMin && tt < tMax)
0172 ener += en;
0173 }
0174 if (debug)
0175 edm::LogVerbatim("IsoTrack") << st1.str();
0176 if (ener > ethr)
0177 energySum += ener;
0178 }
0179 }
0180 if (debug)
0181 edm::LogVerbatim("IsoTrack") << "energyECAL: energySum = " << energySum;
0182 return energySum;
0183 }
0184
0185 template <typename T>
0186 double energyECAL(std::vector<DetId>& vdets,
0187 edm::Handle<T>& hitsEB,
0188 edm::Handle<T>& hitsEE,
0189 const EcalTrigTowerConstituentsMap& ttMap,
0190 double ebThr,
0191 double eeThr,
0192 double tMin,
0193 double tMax,
0194 bool debug) {
0195 double energySum = 0.0;
0196 for (unsigned int i1 = 0; i1 < vdets.size(); i1++) {
0197 if (vdets[i1] != DetId(0)) {
0198 double eTower = spr::energyECALTower(vdets[i1], hitsEB, hitsEE, ttMap, debug);
0199 bool ok = true;
0200 if (vdets[i1].subdetId() == EcalBarrel)
0201 ok = (eTower > ebThr);
0202 else if (vdets[i1].subdetId() == EcalEndcap)
0203 ok = (eTower > eeThr);
0204 std::ostringstream st1;
0205 if (debug)
0206 st1 << "Crystal 0x" << std::hex << vdets[i1]() << std::dec << " Flag " << ok;
0207 if (ok) {
0208 std::vector<typename T::const_iterator> hit;
0209 if (vdets[i1].subdetId() == EcalBarrel) {
0210 spr::findHit(hitsEB, vdets[i1], hit, debug);
0211 } else if (vdets[i1].subdetId() == EcalEndcap) {
0212 spr::findHit(hitsEE, vdets[i1], hit, debug);
0213 }
0214 double ener = 0;
0215 for (unsigned int ihit = 0; ihit < hit.size(); ihit++) {
0216 double en = 0, tt = 0;
0217 if (vdets[i1].subdetId() == EcalBarrel) {
0218 if (hit[ihit] != hitsEB->end()) {
0219 en = hit[ihit]->energy();
0220 tt = hit[ihit]->time();
0221 }
0222 } else if (vdets[i1].subdetId() == EcalEndcap) {
0223 if (hit[ihit] != hitsEE->end()) {
0224 en = hit[ihit]->energy();
0225 tt = hit[ihit]->time();
0226 }
0227 }
0228 if (debug)
0229 st1 << " " << ihit << " E " << en << " time " << tt;
0230 if (tt > tMin && tt < tMax)
0231 ener += en;
0232 }
0233 energySum += ener;
0234 }
0235 if (debug)
0236 edm::LogVerbatim("IsoTrack") << st1.str();
0237 }
0238 }
0239 if (debug)
0240 edm::LogVerbatim("IsoTrack") << "energyECAL: energySum = " << energySum;
0241 return energySum;
0242 }
0243
0244 template <typename T>
0245 double energyECALTower(const DetId& detId,
0246 edm::Handle<T>& hitsEB,
0247 edm::Handle<T>& hitsEE,
0248 const EcalTrigTowerConstituentsMap& ttMap,
0249 bool debug) {
0250 double ener = 0;
0251 EcalTrigTowerDetId trId = ttMap.towerOf(detId);
0252 std::vector<DetId> vdets = ttMap.constituentsOf(trId);
0253 if (debug) {
0254 std::ostringstream st1;
0255 st1 << "energyECALTower: ";
0256 if (detId.subdetId() == EcalBarrel) {
0257 EBDetId id = detId;
0258 st1 << "Cell 0x" << std::hex << detId() << std::dec << " " << id;
0259 } else if (detId.subdetId() == EcalEndcap) {
0260 EEDetId id = detId;
0261 st1 << "Cell 0x" << std::hex << detId() << std::dec << " " << id;
0262 } else {
0263 st1 << "Cell 0x" << std::hex << detId() << std::dec << " Unknown Type";
0264 }
0265 edm::LogVerbatim("IsoTrack") << st1.str() << " Tower " << trId << " with " << vdets.size() << " cells";
0266 }
0267 for (unsigned int i1 = 0; i1 < vdets.size(); i1++) {
0268 if (vdets[i1] != DetId(0)) {
0269 std::vector<typename T::const_iterator> hit;
0270 if (vdets[i1].subdetId() == EcalBarrel) {
0271 spr::findHit(hitsEB, vdets[i1], hit, debug);
0272 } else if (vdets[i1].subdetId() == EcalEndcap) {
0273 spr::findHit(hitsEE, vdets[i1], hit, debug);
0274 }
0275 std::ostringstream st1;
0276 if (debug)
0277 st1 << "Xtal 0x" << std::hex << vdets[i1]() << std::dec;
0278 double en = 0;
0279 for (unsigned int ihit = 0; ihit < hit.size(); ihit++) {
0280 if (vdets[i1].subdetId() == EcalBarrel) {
0281 if (hit[ihit] != hitsEB->end())
0282 en += hit[ihit]->energy();
0283 } else if (vdets[i1].subdetId() == EcalEndcap) {
0284 if (hit[ihit] != hitsEE->end())
0285 en += hit[ihit]->energy();
0286 }
0287 }
0288 if (debug)
0289 edm::LogVerbatim("IsoTrack") << st1.str() << " " << hit.size() << " E " << en;
0290 ener += en;
0291 }
0292 }
0293 if (debug)
0294 edm::LogVerbatim("IsoTrack") << "energyECALTower: Energy in the Tower = " << ener;
0295 return ener;
0296 }
0297
0298 template <typename T>
0299 DetId hotCrystal(const DetId& detId,
0300 edm::Handle<T>& hitsEB,
0301 edm::Handle<T>& hitsEE,
0302 const CaloGeometry* geo,
0303 const CaloTopology* caloTopology,
0304 int ieta,
0305 int iphi,
0306 double tMin,
0307 double tMax,
0308 bool debug) {
0309 std::vector<DetId> vdets;
0310 spr::matrixECALIds(detId, ieta, iphi, geo, caloTopology, vdets, false);
0311
0312 if (debug) {
0313 edm::LogVerbatim("IsoTrack") << "Inside eECALmatrix " << 2 * ieta + 1 << "X" << 2 * iphi + 1 << " nXtals "
0314 << vdets.size();
0315 spr::debugEcalDets(0, vdets);
0316 }
0317
0318 DetId det = spr::hotCrystal(vdets, hitsEB, hitsEE, tMin, tMax, debug);
0319 if (det == DetId(0))
0320 det = detId;
0321 return det;
0322 }
0323
0324 template <typename T>
0325 DetId hotCrystal(
0326 std::vector<DetId>& vdets, edm::Handle<T>& hitsEB, edm::Handle<T>& hitsEE, double tMin, double tMax, bool debug) {
0327 DetId det = DetId(0);
0328 double eMax = -99999.;
0329 for (unsigned int i1 = 0; i1 < vdets.size(); i1++) {
0330 if (vdets[i1] != DetId(0)) {
0331 std::vector<typename T::const_iterator> hit;
0332 if (vdets[i1].subdetId() == EcalBarrel) {
0333 spr::findHit(hitsEB, vdets[i1], hit, debug);
0334 } else if (vdets[i1].subdetId() == EcalEndcap) {
0335 spr::findHit(hitsEE, vdets[i1], hit, debug);
0336 }
0337 std::ostringstream st1;
0338 if (debug)
0339 st1 << "Crystal 0x" << std::hex << vdets[i1]() << std::dec;
0340 double ener = 0;
0341 for (unsigned int ihit = 0; ihit < hit.size(); ihit++) {
0342 double en = 0, tt = 0;
0343 if (vdets[i1].subdetId() == EcalBarrel) {
0344 if (hit[ihit] != hitsEB->end()) {
0345 en = hit[ihit]->energy();
0346 tt = hit[ihit]->time();
0347 }
0348 } else if (vdets[i1].subdetId() == EcalEndcap) {
0349 if (hit[ihit] != hitsEE->end()) {
0350 en = hit[ihit]->energy();
0351 tt = hit[ihit]->time();
0352 }
0353 }
0354 if (debug)
0355 st1 << " " << ihit << " " << en << " " << tt;
0356 if (tt > tMin && tt < tMax)
0357 ener += en;
0358 }
0359 if (debug)
0360 edm::LogVerbatim("IsoTrack") << st1.str();
0361 if (ener > eMax) {
0362 det = vdets[i1];
0363 eMax = ener;
0364 }
0365 }
0366 }
0367 if (debug)
0368 edm::LogVerbatim("IsoTrack") << "energyECAL: maxEnegy = " << eMax;
0369 return det;
0370 }
0371
0372 template <typename T>
0373 double eECALmatrix(CaloNavigator<DetId>& navigator,
0374 edm::Handle<T>& hits,
0375 int ieta,
0376 int iphi,
0377 const EcalSeverityLevelAlgo* sevlv,
0378 bool debug) {
0379 std::vector<typename T::const_iterator> hit;
0380 spr::hitECALmatrix(navigator, hits, ieta, iphi, hit, debug);
0381
0382 if (debug) {
0383 edm::LogVerbatim("IsoTrack") << "Inside eECALmatrix " << 2 * ieta + 1 << "X" << 2 * iphi + 1;
0384 std::ostringstream st1;
0385 st1 << "nXtals " << hit.size();
0386 for (unsigned int ihit = 0; ihit < hit.size(); ihit++)
0387 st1 << " ihit:" << ihit << " " << (unsigned int)hit[ihit]->id();
0388 edm::LogVerbatim("IsoTrack") << st1.str();
0389 }
0390
0391 double energySum = 0.0;
0392 for (unsigned int ihit = 0; ihit < hit.size(); ihit++) {
0393 energySum += hit[ihit]->energy();
0394 }
0395 return energySum;
0396 }
0397
0398 template <typename T>
0399 std::vector<std::pair<DetId, double> > eECALmatrixCell(const DetId& detId,
0400 edm::Handle<T>& hitsEB,
0401 edm::Handle<T>& hitsEE,
0402 const CaloGeometry* geo,
0403 const CaloTopology* caloTopology,
0404 int ieta,
0405 int iphi,
0406 double ebThr,
0407 double eeThr,
0408 bool debug) {
0409 std::vector<DetId> vdets = spr::matrixECALIds(detId, ieta, iphi, geo, caloTopology, debug);
0410
0411 if (debug) {
0412 edm::LogVerbatim("IsoTrack") << "Inside eECALmatrixCell " << 2 * ieta + 1 << "X" << 2 * iphi + 1 << " nXtals "
0413 << vdets.size();
0414 }
0415
0416 return spr::energyECALCell(vdets, hitsEB, hitsEE, ebThr, eeThr, debug);
0417 }
0418
0419 template <typename T>
0420 std::pair<double, int> eECALmatrixTotal(const DetId& detId,
0421 edm::Handle<T>& hitsEB,
0422 edm::Handle<T>& hitsEE,
0423 const CaloGeometry* geo,
0424 const CaloTopology* caloTopology,
0425 int ieta,
0426 int iphi,
0427 double ebThr,
0428 double eeThr,
0429 double tMin,
0430 double tMax,
0431 bool debug) {
0432 spr::EtaPhi etaphi = spr::getEtaPhi(ieta, iphi, debug);
0433
0434
0435 double energySum = 0;
0436 int itrym = 0;
0437 for (int itry = 0; itry < etaphi.ntrys; itry++) {
0438 std::vector<DetId> vdets = spr::matrixECALIds(detId,
0439 etaphi.ietaE[itry],
0440 etaphi.ietaW[itry],
0441 etaphi.iphiN[itry],
0442 etaphi.iphiS[itry],
0443 geo,
0444 caloTopology,
0445 debug);
0446 double energy = spr::energyECAL(vdets, hitsEB, hitsEE, ebThr, eeThr, tMin, tMax, debug);
0447 if (energy > energySum) {
0448 energySum = energy;
0449 itrym = itry;
0450 }
0451 }
0452
0453 if (debug)
0454 edm::LogVerbatim("IsoTrack") << "eECALmatrixTotal:: energy deposit in " << ieta << "X" << iphi << " matrix is "
0455 << energySum << " for trial # " << itrym;
0456 return std::pair<double, int>(energySum, itrym);
0457 }
0458
0459 template <typename T>
0460 std::vector<std::pair<DetId, double> > energyECALCell(std::vector<DetId>& vdets,
0461 edm::Handle<T>& hitsEB,
0462 edm::Handle<T>& hitsEE,
0463 double ebThr,
0464 double eeThr,
0465 bool debug) {
0466 std::vector<std::pair<DetId, double> > energyDets;
0467 for (unsigned int i1 = 0; i1 < vdets.size(); i1++) {
0468 double energySum = 0.0;
0469 if (vdets[i1] != DetId(0)) {
0470 std::vector<typename T::const_iterator> hit;
0471 if (vdets[i1].subdetId() == EcalBarrel) {
0472 hit = spr::findHit(hitsEB, vdets[i1]);
0473 } else if (vdets[i1].subdetId() == EcalEndcap) {
0474 hit = spr::findHit(hitsEE, vdets[i1]);
0475 }
0476 std::ostringstream st1;
0477 if (debug)
0478 st1 << "Xtal 0x" << std::hex << vdets[i1]() << std::dec;
0479 bool ok = false;
0480 double ethr = ebThr;
0481 for (unsigned int ihit = 0; ihit < hit.size(); ihit++) {
0482 double en = 0;
0483 if (vdets[i1].subdetId() == EcalBarrel) {
0484 if (hit[ihit] != hitsEB->end()) {
0485 en = hit[ihit]->energy();
0486 ok = true;
0487 }
0488 } else if (vdets[i1].subdetId() == EcalEndcap) {
0489 if (hit[ihit] != hitsEE->end()) {
0490 en = hit[ihit]->energy();
0491 ok = true;
0492 ethr = eeThr;
0493 }
0494 }
0495 if (debug)
0496 st1 << " " << ihit << " " << en;
0497 energySum += en;
0498 }
0499 if (debug)
0500 edm::LogVerbatim("IsoTrack") << st1.str() << "\nenergyECALCell: energySum = " << energySum;
0501 if (ok && energySum > ethr)
0502 energyDets.push_back(std::pair<DetId, double>(vdets[i1], energySum));
0503 }
0504 }
0505 return energyDets;
0506 }
0507
0508 }
