Back to home page

Project CMSSW displayed by LXR

 
 

    


File indexing completed on 2024-04-06 11:59:20

0001 #include "Calibration/IsolatedParticles/interface/eECALMatrix.h"
0002 #include "Calibration/IsolatedParticles/interface/FindCaloHit.h"
0003 #include "DataFormats/CaloTowers/interface/CaloTowerDetId.h"
0004 #include "DataFormats/EcalDetId/interface/EBDetId.h"
0005 #include "DataFormats/EcalDetId/interface/EEDetId.h"
0006 #include "DataFormats/EcalDetId/interface/ESDetId.h"
0007 #include "DataFormats/EcalDetId/interface/EcalSubdetector.h"
0008 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0009 #include "RecoLocalCalo/EcalRecAlgos/interface/EcalSeverityLevelAlgo.h"
0010 
0011 #include <sstream>
0012 
0013 namespace spr {
0014 
0015   std::pair<double, bool> energyECAL(const DetId& id,
0016                                      edm::Handle<EcalRecHitCollection>& hitsEC,
0017                                      const EcalSeverityLevelAlgo* sevlv,
0018                                      bool testSpike,
0019                                      double tMin,
0020                                      double tMax,
0021                                      bool debug) {
0022     std::vector<EcalRecHitCollection::const_iterator> hits;
0023     spr::findHit(hitsEC, id, hits, debug);
0024 
0025     std::ostringstream st1;
0026     if (debug)
0027       st1 << "Xtal 0x" << std::hex << id() << std::dec;
0028     const EcalRecHitCollection* recHitsEC = (hitsEC.isValid()) ? hitsEC.product() : nullptr;
0029     bool flag = (!testSpike) ? true : (sevlv->severityLevel(id, (*recHitsEC)) != EcalSeverityLevel::kWeird);
0030     double ener(0);
0031     for (const auto& hit : hits) {
0032       double en(0), tt(0);
0033       if (hit != hitsEC->end()) {
0034         en = hit->energy();
0035         tt = hit->time();
0036       }
0037       if (debug)
0038         st1 << " " << tt << " " << en;
0039       if (tt > tMin && tt < tMax)
0040         ener += en;
0041     }
0042     if (debug) {
0043       if (!flag)
0044         st1 << " detected to be a spike";
0045       edm::LogVerbatim("IsoTrack") << st1.str();
0046     }
0047     return std::pair<double, bool>(ener, flag);
0048   }
0049 
0050   std::pair<double, bool> energyECAL(const std::vector<DetId>& vdets,
0051                                      edm::Handle<EcalRecHitCollection>& hitsEC,
0052                                      const EcalSeverityLevelAlgo* sevlv,
0053                                      bool noThrCut,
0054                                      bool testSpike,
0055                                      double eThr,
0056                                      double tMin,
0057                                      double tMax,
0058                                      bool debug) {
0059     bool flag(true);
0060     double energySum(0.0);
0061     for (const auto& id : vdets) {
0062       if (id != DetId(0)) {
0063         std::pair<double, bool> ecalEn = spr::energyECAL(id, hitsEC, sevlv, testSpike, tMin, tMax, debug);
0064         if (!ecalEn.second)
0065           flag = false;
0066         if ((ecalEn.first > eThr) || noThrCut)
0067           energySum += ecalEn.first;
0068       }
0069     }
0070     if (debug)
0071       edm::LogVerbatim("IsoTrack") << "energyECAL: energySum = " << energySum << " flag = " << flag;
0072     return std::pair<double, bool>(energySum, flag);
0073   }
0074 
0075   std::pair<double, bool> eECALmatrix(const DetId& detId,
0076                                       edm::Handle<EcalRecHitCollection>& hitsEB,
0077                                       edm::Handle<EcalRecHitCollection>& hitsEE,
0078                                       const EcalChannelStatus& chStatus,
0079                                       const CaloGeometry* geo,
0080                                       const CaloTopology* caloTopology,
0081                                       const EcalSeverityLevelAlgo* sevlv,
0082                                       int ieta,
0083                                       int iphi,
0084                                       double ebThr,
0085                                       double eeThr,
0086                                       double tMin,
0087                                       double tMax,
0088                                       bool debug) {
0089     std::vector<DetId> vdets;
0090     spr::matrixECALIds(detId, ieta, iphi, geo, caloTopology, vdets, debug);
0091     if (debug)
0092       edm::LogVerbatim("IsoTrack") << "Inside eECALmatrix " << 2 * ieta + 1 << "X" << 2 * iphi + 1 << " nXtals "
0093                                    << vdets.size();
0094     bool flag(true);
0095     for (const auto& id : vdets) {
0096       if ((id.det() == DetId::Ecal) && (id.subdetId() == EcalBarrel)) {
0097         if (sevlv->severityLevel(id, (*hitsEB)) == EcalSeverityLevel::kWeird)
0098           flag = false;
0099       } else if ((id.det() == DetId::Ecal) && (id.subdetId() == EcalEndcap)) {
0100         if (sevlv->severityLevel(id, (*hitsEE)) == EcalSeverityLevel::kWeird)
0101           flag = false;
0102       }
0103     }
0104     return std::pair<double, bool>(spr::energyECAL(vdets, hitsEB, hitsEE, ebThr, eeThr, tMin, tMax, debug), flag);
0105   }
0106 
0107   std::pair<double, bool> eECALmatrix(const DetId& detId,
0108                                       edm::Handle<EcalRecHitCollection>& hitsEB,
0109                                       edm::Handle<EcalRecHitCollection>& hitsEE,
0110                                       const EcalChannelStatus& chStatus,
0111                                       const CaloGeometry* geo,
0112                                       const CaloTopology* caloTopology,
0113                                       const EcalSeverityLevelAlgo* sevlv,
0114                                       const EcalTrigTowerConstituentsMap& ttMap,
0115                                       int ieta,
0116                                       int iphi,
0117                                       double ebThr,
0118                                       double eeThr,
0119                                       double tMin,
0120                                       double tMax,
0121                                       bool debug) {
0122     std::vector<DetId> vdets;
0123     spr::matrixECALIds(detId, ieta, iphi, geo, caloTopology, vdets, debug);
0124     if (debug)
0125       edm::LogVerbatim("IsoTrack") << "Inside eECALmatrix " << 2 * ieta + 1 << "X" << 2 * iphi + 1 << " nXtals "
0126                                    << vdets.size();
0127 
0128     bool flag(true);
0129     double energySum = 0.0;
0130     for (const auto& id : vdets) {
0131       if ((id != DetId(0)) && (id.det() == DetId::Ecal) &&
0132           ((id.subdetId() == EcalBarrel) || (id.subdetId() == EcalEndcap))) {
0133         double eTower = spr::energyECALTower(id, hitsEB, hitsEE, ttMap, debug);
0134         bool ok = (id.subdetId() == EcalBarrel) ? (eTower > ebThr) : (eTower > eeThr);
0135         if (debug && (!ok))
0136           edm::LogVerbatim("IsoTrack") << "Crystal 0x" << std::hex << id() << std::dec << " Flag " << ok;
0137         if (ok) {
0138           std::pair<double, bool> ecalEn = (id.subdetId() == EcalBarrel)
0139                                                ? spr::energyECAL(id, hitsEB, sevlv, true, tMin, tMax, debug)
0140                                                : spr::energyECAL(id, hitsEE, sevlv, false, tMin, tMax, debug);
0141           if (!ecalEn.second)
0142             flag = false;
0143           energySum += ecalEn.first;
0144         }
0145       }
0146     }
0147     if (debug)
0148       edm::LogVerbatim("IsoTrack") << "energyECAL: energySum = " << energySum << " flag = " << flag;
0149     return std::pair<double, bool>(energySum, flag);
0150   }
0151 
0152   std::pair<double, bool> eECALmatrix(const HcalDetId& detId,
0153                                       edm::Handle<EcalRecHitCollection>& hitsEB,
0154                                       edm::Handle<EcalRecHitCollection>& hitsEE,
0155                                       const CaloGeometry* geo,
0156                                       const CaloTowerConstituentsMap* ctmap,
0157                                       const EcalSeverityLevelAlgo* sevlv,
0158                                       double ebThr,
0159                                       double eeThr,
0160                                       double tMin,
0161                                       double tMax,
0162                                       bool debug) {
0163     CaloTowerDetId tower = ctmap->towerOf(detId);
0164     std::vector<DetId> ids = ctmap->constituentsOf(tower);
0165     if (debug) {
0166       edm::LogVerbatim("IsoTrack") << "eECALmatrix: " << detId << " belongs to " << tower << " which has " << ids.size()
0167                                    << " constituents";
0168       for (unsigned int i = 0; i < ids.size(); ++i) {
0169         std::ostringstream st1;
0170         st1 << "[" << i << "] " << std::hex << ids[i].rawId() << std::dec;
0171         if (ids[i].det() == DetId::Ecal && ids[i].subdetId() == EcalBarrel) {
0172           st1 << " " << EBDetId(ids[i]);
0173         } else if (ids[i].det() == DetId::Ecal && ids[i].subdetId() == EcalEndcap) {
0174           st1 << " " << EEDetId(ids[i]);
0175         } else if (ids[i].det() == DetId::Ecal && ids[i].subdetId() == EcalPreshower) {
0176           st1 << " " << ESDetId(ids[i]);
0177         } else if (ids[i].det() == DetId::Hcal) {
0178           st1 << " " << HcalDetId(ids[i]);
0179         }
0180         edm::LogVerbatim("IsoTrack") << st1.str();
0181       }
0182     }
0183 
0184     std::vector<DetId> idEBEE;
0185     bool flag(true);
0186     for (const auto& id : ids) {
0187       if ((id.det() == DetId::Ecal) && (id.subdetId() == EcalBarrel)) {
0188         idEBEE.emplace_back(id);
0189         if (sevlv->severityLevel(id, (*hitsEB)) == EcalSeverityLevel::kWeird)
0190           flag = false;
0191       } else if ((id.det() == DetId::Ecal) && (id.subdetId() == EcalEndcap)) {
0192         idEBEE.emplace_back(id);
0193         if (sevlv->severityLevel(id, (*hitsEE)) == EcalSeverityLevel::kWeird)
0194           flag = false;
0195       }
0196     }
0197 
0198     if (debug)
0199       edm::LogVerbatim("IsoTrack") << "eECALmatrix: with " << idEBEE.size() << " EB+EE hits and "
0200                                    << "spike flag " << flag;
0201     double etot = (!idEBEE.empty()) ? spr::energyECAL(idEBEE, hitsEB, hitsEE, ebThr, eeThr, tMin, tMax, debug) : 0;
0202     return std::pair<double, bool>(etot, flag);
0203   }
0204 
0205   std::pair<double, bool> eECALmatrix(const DetId& detId,
0206                                       edm::Handle<EcalRecHitCollection>& hitsEB,
0207                                       edm::Handle<EcalRecHitCollection>& hitsEE,
0208                                       const EcalChannelStatus& chStatus,
0209                                       const CaloGeometry* geo,
0210                                       const CaloTopology* caloTopology,
0211                                       const EcalSeverityLevelAlgo* sevlv,
0212                                       const EcalPFRecHitThresholds* eThresholds,
0213                                       int ieta,
0214                                       int iphi,
0215                                       bool debug) {
0216     std::vector<DetId> vdets;
0217     spr::matrixECALIds(detId, ieta, iphi, geo, caloTopology, vdets, debug);
0218     if (debug)
0219       edm::LogVerbatim("IsoTrack") << "Inside eECALmatrix " << 2 * ieta + 1 << "X" << 2 * iphi + 1 << " nXtals "
0220                                    << vdets.size();
0221     bool flag(true);
0222     for (const auto& id : vdets) {
0223       if ((id.det() == DetId::Ecal) && (id.subdetId() == EcalBarrel)) {
0224         if (sevlv->severityLevel(id, (*hitsEB)) == EcalSeverityLevel::kWeird)
0225           flag = false;
0226       } else if ((id.det() == DetId::Ecal) && (id.subdetId() == EcalEndcap)) {
0227         if (sevlv->severityLevel(id, (*hitsEE)) == EcalSeverityLevel::kWeird)
0228           flag = false;
0229       }
0230     }
0231     double energySum = 0.0;
0232     for (unsigned int i1 = 0; i1 < vdets.size(); i1++) {
0233       if (vdets[i1] != DetId(0)) {
0234         std::vector<EcalRecHitCollection::const_iterator> hit;
0235         if (vdets[i1].subdetId() == EcalBarrel) {
0236           spr::findHit(hitsEB, vdets[i1], hit, debug);
0237         } else if (vdets[i1].subdetId() == EcalEndcap) {
0238           spr::findHit(hitsEE, vdets[i1], hit, debug);
0239         }
0240         std::ostringstream st1;
0241         if (debug)
0242           st1 << "Crystal 0x" << std::hex << vdets[i1]() << std::dec;
0243         double ener = 0, ethr = static_cast<double>((*eThresholds)[vdets[i1]]);
0244         for (unsigned int ihit = 0; ihit < hit.size(); ihit++) {
0245           double en = 0;
0246           if (vdets[i1].subdetId() == EcalBarrel) {
0247             if (hit[ihit] != hitsEB->end())
0248               en = hit[ihit]->energy();
0249           } else if (vdets[i1].subdetId() == EcalEndcap) {
0250             if (hit[ihit] != hitsEE->end())
0251               en = hit[ihit]->energy();
0252           }
0253           if (debug)
0254             st1 << " " << ihit << " " << en << " Thr " << ethr;
0255           ener += en;
0256         }
0257         if (debug)
0258           edm::LogVerbatim("IsoTrack") << st1.str();
0259         if (ener > ethr)
0260           energySum += ener;
0261       }
0262     }
0263     return std::pair<double, bool>(energySum, flag);
0264   }
0265 
0266 }  // namespace spr