Back to home page

Project CMSSW displayed by LXR

 
 

    


File indexing completed on 2024-09-07 04:35:02

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         // shift the navigator by dx/dy crystals in eta/phi
0097         thisDet = navigator.offsetBy(dx, dy);
0098 
0099         // Place the navigator back to the original position
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       }  // iphi
0129     }  // ieta
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     // Get maximum of all the trials
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 }  // namespace spr