Back to home page

Project CMSSW displayed by LXR

 
 

    


File indexing completed on 2021-08-12 03:21:13

0001 #include "Calibration/IsolatedParticles/interface/DebugInfo.h"
0002 #include "Calibration/IsolatedParticles/interface/eHCALMatrix.h"
0003 #include "DataFormats/HcalDetId/interface/HcalSubdetector.h"
0004 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0005 
0006 #include <algorithm>
0007 #include <sstream>
0008 
0009 namespace spr {
0010   double eHCALmatrix(const HcalTopology* topology,
0011                      const DetId& det0,
0012                      std::vector<PCaloHit>& hits,
0013                      int ieta,
0014                      int iphi,
0015                      bool includeHO,
0016                      double hbThr,
0017                      double heThr,
0018                      double hfThr,
0019                      double hoThr,
0020                      double tMin,
0021                      double tMax,
0022                      bool debug) {
0023     HcalDetId hcid0(det0.rawId());
0024     HcalDetId hcid(hcid0.subdet(), hcid0.ieta(), hcid0.iphi(), 1);
0025     DetId det(hcid.rawId());
0026     if (debug)
0027       edm::LogVerbatim("IsoTrack") << "Inside eHCALmatrix " << 2 * ieta + 1 << "X" << 2 * iphi + 1
0028                                    << " Inclusion of HO Flag " << includeHO;
0029 
0030     double energySum(0);
0031     std::vector<DetId> dets(1, det);
0032     std::vector<DetId> vdets = spr::matrixHCALIds(dets, topology, ieta, iphi, includeHO, false);
0033     if (debug) {
0034       edm::LogVerbatim("IsoTrack") << "matrixHCALIds::Total number of cells found is " << vdets.size();
0035       spr::debugHcalDets(0, vdets);
0036     }
0037 
0038     int khit(0);
0039     for (unsigned int i = 0; i < vdets.size(); i++) {
0040       std::vector<std::vector<PCaloHit>::const_iterator> hit = spr::findHit(hits, vdets[i]);
0041       double energy = 0;
0042       int subdet = ((HcalDetId)(vdets[i].rawId())).subdet();
0043       double eThr = spr::eHCALThreshold(subdet, hbThr, heThr, hfThr, hoThr);
0044       for (unsigned int ihit = 0; ihit < hit.size(); ihit++) {
0045         if (hit[ihit] != hits.end()) {
0046           khit++;
0047           if (debug)
0048             edm::LogVerbatim("IsoTrack") << "energyHCAL:: Hit " << khit << " " << (HcalDetId)vdets[i] << " E "
0049                                          << hit[ihit]->energy() << " t " << hit[ihit]->time();
0050 
0051           if (hit[ihit]->time() > tMin && hit[ihit]->time() < tMax) {
0052             energy += hit[ihit]->energy();
0053           }
0054         }
0055       }
0056       if (energy > eThr)
0057         energySum += energy;
0058     }
0059 
0060     if (debug)
0061       edm::LogVerbatim("IsoTrack") << "eHCALmatrix::Total energy " << energySum;
0062     return energySum;
0063   }
0064 
0065   double eHCALmatrix(const CaloGeometry* geo,
0066                      const HcalTopology* topology,
0067                      const DetId& det0,
0068                      std::vector<PCaloHit>& hits,
0069                      int ieta,
0070                      int iphi,
0071                      HcalDetId& hotCell,
0072                      bool includeHO,
0073                      bool debug) {
0074     HcalDetId hcid0(det0.rawId());
0075     HcalDetId hcid(hcid0.subdet(), hcid0.ieta(), hcid0.iphi(), 1);
0076     DetId det(hcid.rawId());
0077     std::vector<DetId> dets(1, det);
0078     std::vector<DetId> vdets = spr::matrixHCALIds(dets, topology, ieta, iphi, includeHO, debug);
0079     hotCell = hcid0;
0080 
0081     std::vector<std::vector<PCaloHit>::const_iterator> hitlist;
0082     for (unsigned int i = 0; i < vdets.size(); i++) {
0083       std::vector<std::vector<PCaloHit>::const_iterator> hit = spr::findHit(hits, vdets[i]);
0084       hitlist.insert(hitlist.end(), hit.begin(), hit.end());
0085     }
0086 
0087     double energySum(0);
0088     for (unsigned int ihit = 0; ihit < hitlist.size(); ihit++)
0089       energySum += hitlist[ihit]->energy();
0090 
0091     // Get hotCell ID
0092     dets.clear();
0093     std::vector<double> energies;
0094     for (unsigned int ihit = 0; ihit < hitlist.size(); ihit++) {
0095       double energy = hitlist[ihit]->energy();
0096       HcalDetId id0 = HcalDetId(hitlist[ihit]->id());
0097       if ((id0.subdet() != HcalOuter) || includeHO) {
0098         HcalDetId id1(id0.subdet(), id0.ieta(), id0.iphi(), 1);
0099         bool found(false);
0100         for (unsigned int idet = 0; idet < dets.size(); ++idet) {
0101           if (id1 == HcalDetId(dets[idet])) {
0102             energies[idet] += energy;
0103             found = true;
0104             break;
0105           }
0106         }
0107         if (!found) {
0108           dets.push_back(DetId(id1));
0109           energies.push_back(energy);
0110         }
0111       }
0112     }
0113     double energyMax(-99.);
0114     for (unsigned int ihit = 0; ihit < dets.size(); ihit++) {
0115       if (energies[ihit] > energyMax) {
0116         energyMax = energies[ihit];
0117         hotCell = HcalDetId(dets[ihit]);
0118       }
0119     }
0120     return energySum;
0121   }
0122 
0123   void energyHCALCell(HcalDetId detId,
0124                       std::vector<PCaloHit>& hits,
0125                       std::vector<std::pair<double, int> >& energyCell,
0126                       int maxDepth,
0127                       double hbThr,
0128                       double heThr,
0129                       double hfThr,
0130                       double hoThr,
0131                       double tMin,
0132                       double tMax,
0133                       int depthHE,
0134                       bool debug) {
0135     energyCell.clear();
0136     int subdet = detId.subdet();
0137     double eThr = spr::eHCALThreshold(subdet, hbThr, heThr, hfThr, hoThr);
0138     bool hbhe = (detId.ietaAbs() == 16);
0139     if (debug)
0140       edm::LogVerbatim("IsoTrack") << "energyHCALCell: input ID " << detId << " MaxDepth " << maxDepth
0141                                    << " Threshold (E) " << eThr << " (T) " << tMin << ":" << tMax;
0142 
0143     for (int i = 0; i < maxDepth; i++) {
0144       HcalSubdetector subdet0 = (hbhe) ? ((i + 1 >= depthHE) ? HcalEndcap : HcalBarrel) : detId.subdet();
0145       HcalDetId hcid(subdet0, detId.ieta(), detId.iphi(), i + 1);
0146       DetId det(hcid.rawId());
0147       std::vector<std::vector<PCaloHit>::const_iterator> hit = spr::findHit(hits, det);
0148       double energy(0);
0149       for (unsigned int ihit = 0; ihit < hit.size(); ++ihit) {
0150         if (hit[ihit]->time() > tMin && hit[ihit]->time() < tMax)
0151           energy += hit[ihit]->energy();
0152         if (debug)
0153           edm::LogVerbatim("IsoTrack") << "energyHCALCell:: Hit[" << ihit << "] " << hcid << " E "
0154                                        << hit[ihit]->energy() << " t " << hit[ihit]->time();
0155       }
0156 
0157       if (debug)
0158         edm::LogVerbatim("IsoTrack") << "energyHCALCell:: Cell " << hcid << " E " << energy << " from " << hit.size()
0159                                      << " threshold " << eThr;
0160       if (energy > eThr && !hit.empty()) {
0161         energyCell.push_back(std::pair<double, int>(energy, i + 1));
0162       }
0163     }
0164 
0165     if (debug) {
0166       std::ostringstream st1;
0167       st1 << "energyHCALCell:: " << energyCell.size() << " entries from " << maxDepth << " depths:";
0168       for (unsigned int i = 0; i < energyCell.size(); ++i) {
0169         st1 << " [" << i << "] (" << energyCell[i].first << ":" << energyCell[i].second << ")";
0170       }
0171       edm::LogVerbatim("IsoTrack") << st1.str();
0172     }
0173   }
0174 
0175   HcalDetId getHotCell(std::vector<HBHERecHitCollection::const_iterator>& hit, bool includeHO, int useRaw, bool) {
0176     std::vector<HcalDetId> dets;
0177     std::vector<double> energies;
0178     for (unsigned int ihit = 0; ihit < hit.size(); ihit++) {
0179       double energy = getRawEnergy(hit.at(ihit), useRaw);
0180       HcalDetId id0 = hit.at(ihit)->id();
0181       if ((id0.subdet() != HcalOuter) || includeHO) {
0182         HcalDetId id1(id0.subdet(), id0.ieta(), id0.iphi(), 1);
0183         bool found(false);
0184         for (unsigned int idet = 0; idet < dets.size(); ++idet) {
0185           if (id1 == dets[idet]) {
0186             energies[idet] += energy;
0187             found = true;
0188             break;
0189           }
0190         }
0191         if (!found) {
0192           dets.push_back(id1);
0193           energies.push_back(energy);
0194         }
0195       }
0196     }
0197     double energyMax(-99.);
0198     HcalDetId hotCell;
0199     for (unsigned int ihit = 0; ihit < dets.size(); ihit++) {
0200       if (energies[ihit] > energyMax) {
0201         energyMax = energies[ihit];
0202         hotCell = dets[ihit];
0203       }
0204     }
0205     return hotCell;
0206   }
0207 
0208   HcalDetId getHotCell(std::vector<std::vector<PCaloHit>::const_iterator>& hit, bool includeHO, int useRaw, bool) {
0209     std::vector<HcalDetId> dets;
0210     std::vector<double> energies;
0211     for (unsigned int ihit = 0; ihit < hit.size(); ihit++) {
0212       double energy = hit.at(ihit)->energy();
0213       HcalDetId id0 = getRawEnergy(hit.at(ihit), useRaw);
0214       if ((id0.subdet() != HcalOuter) || includeHO) {
0215         HcalDetId id1(id0.subdet(), id0.ieta(), id0.iphi(), 1);
0216         bool found(false);
0217         for (unsigned int idet = 0; idet < dets.size(); ++idet) {
0218           if (id1 == dets[idet]) {
0219             energies[idet] += energy;
0220             found = true;
0221             break;
0222           }
0223         }
0224         if (!found) {
0225           dets.push_back(id1);
0226           energies.push_back(energy);
0227         }
0228       }
0229     }
0230     double energyMax(-99.);
0231     HcalDetId hotCell;
0232     for (unsigned int ihit = 0; ihit < dets.size(); ihit++) {
0233       if (energies[ihit] > energyMax) {
0234         energyMax = energies[ihit];
0235         hotCell = dets[ihit];
0236       }
0237     }
0238     return hotCell;
0239   }
0240 
0241   double eHCALThreshold(int subdet, double hbThr, double heThr, double hfThr, double hoThr) {
0242     double eThr = hbThr;
0243     if (subdet == (int)(HcalEndcap))
0244       eThr = heThr;
0245     else if (subdet == (int)(HcalForward))
0246       eThr = hfThr;
0247     else if (subdet == (int)(HcalOuter))
0248       eThr = hoThr;
0249     return eThr;
0250   }
0251 }  // namespace spr