Back to home page

Project CMSSW displayed by LXR
 
 

    

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 "Calibration/IsolatedParticles/interface/FindDistCone.h"
0005 #include "Calibration/IsolatedParticles/interface/MatrixHCALDetIds.h"
0006 #include "DataFormats/HcalDetId/interface/HcalSubdetector.h"
0007 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0008 
0009 #include <algorithm>
0010 #include <iostream>
0011 #include <sstream>
0012 
0013 namespace spr {
0014 
0015   //*** Returns energy in NxN Hcal towers
0016   template <typename T>
0017   double eHCALmatrix(const HcalTopology* topology,
0018                      const DetId& det0,
0019                      edm::Handle<T>& hits,
0020                      int ieta,
0021                      int iphi,
0022                      bool includeHO,
0023                      bool algoNew,
0024                      double hbThr,
0025                      double heThr,
0026                      double hfThr,
0027                      double hoThr,
0028                      double tMin,
0029                      double tMax,
0030                      int useRaw,
0031                      bool debug) {
0032     std::vector<typename T::const_iterator> hit;
0033     HcalDetId hcid0(det0.rawId());
0034     HcalDetId hcid(hcid0.subdet(), hcid0.ieta(), hcid0.iphi(), 1);
0035     DetId det(hcid.rawId());
0036     if (debug)
0037       edm::LogVerbatim("IsoTrack") << "Inside eHCALmatrix " << 2 * ieta + 1 << "X" << 2 * iphi + 1 << " AlgoNew "
0038                                    << algoNew << " Inclusion of HO Flag " << includeHO;
0039     double energySum = 0.0;
0040     if (algoNew) {
0041       energySum = spr::energyHCALmatrixNew(
0042           topology, det, hits, ieta, iphi, includeHO, hbThr, heThr, hfThr, hoThr, tMin, tMax, useRaw, debug);
0043     } else {
0044       edm::LogVerbatim("IsoTrack") << "eHCALmatrix:: Algorithm New = " << algoNew
0045                                    << " not supported - call directly spr::energyHCALmatrix ";
0046       //       energySum = spr::energyHCALmatrix(topology, det, hits, ieta, iphi, includeHO, hbThr, heThr, hfThr, hoThr, tMin, tMax, useRaw, debug);
0047     }
0048 
0049     if (debug)
0050       edm::LogVerbatim("IsoTrack") << "eHCALmatrix::Total energy " << energySum;
0051     return energySum;
0052   }
0053 
0054   template <typename T>
0055   double eHCALmatrix(const HcalTopology* topology,
0056                      const DetId& det0,
0057                      edm::Handle<T>& hits,
0058                      int ietaE,
0059                      int ietaW,
0060                      int iphiN,
0061                      int iphiS,
0062                      bool includeHO,
0063                      double hbThr,
0064                      double heThr,
0065                      double hfThr,
0066                      double hoThr,
0067                      double tMin,
0068                      double tMax,
0069                      int useRaw,
0070                      bool debug) {
0071     HcalDetId hcid0(det0.rawId());
0072     HcalDetId hcid(hcid0.subdet(), hcid0.ieta(), hcid0.iphi(), 1);
0073     DetId det(hcid.rawId());
0074     if (debug)
0075       edm::LogVerbatim("IsoTrack") << "Inside eHCALmatrix " << ietaE + ietaW + 1 << "X" << iphiN + iphiS + 1
0076                                    << " Inclusion of HO Flag " << includeHO;
0077 
0078     return spr::energyHCALmatrixTotal(topology,
0079                                       det,
0080                                       hits,
0081                                       ietaE,
0082                                       ietaW,
0083                                       iphiN,
0084                                       iphiS,
0085                                       includeHO,
0086                                       hbThr,
0087                                       heThr,
0088                                       hfThr,
0089                                       hoThr,
0090                                       tMin,
0091                                       tMax,
0092                                       useRaw,
0093                                       debug);
0094   }
0095 
0096   // NxN method modified to return rechit vectors for hitmaps
0097   template <typename T>
0098   double eHCALmatrix(const HcalTopology* topology,
0099                      const DetId& det0,
0100                      edm::Handle<T>& hits,
0101                      int ieta,
0102                      int iphi,
0103                      int& nRecHits,
0104                      std::vector<int>& RH_ieta,
0105                      std::vector<int>& RH_iphi,
0106                      std::vector<double>& RH_ene,
0107                      std::set<int>& uniqueIdset,
0108                      int useRaw) {
0109     bool debug = false;
0110     bool includeHO = false;
0111     HcalDetId hcid0(det0.rawId());
0112     HcalDetId hcid(hcid0.subdet(), hcid0.ieta(), hcid0.iphi(), 1);
0113     DetId det(hcid.rawId());
0114 
0115     std::vector<typename T::const_iterator> hit;
0116     spr::hitHCALmatrix(topology, det, hits, ieta, iphi, hit, includeHO, debug);
0117 
0118     nRecHits = -99;
0119     double energySum = 0.0;
0120     for (unsigned int ihit = 0; ihit < hit.size(); ihit++) {
0121       energySum += spr::getEnergy(hit[ihit], useRaw);
0122       int eta(-99);
0123       int phi(-99);
0124       // Get eta phi for counting unique cells
0125       spr::getEtaPhi(hit[ihit], eta, phi);
0126       // Put eta phi in vectors
0127       RH_ieta.push_back(hit[ihit]->id().ieta());
0128       RH_iphi.push_back(hit[ihit]->id().iphi());
0129       RH_ene.push_back(hit[ihit]->energy());
0130       int uniqueEtaPhiId = 100 * eta + phi;
0131       uniqueIdset.insert(uniqueEtaPhiId);
0132     }
0133     nRecHits = uniqueIdset.size();
0134     return energySum;
0135   }
0136 
0137   // NxN method modified to return rechit vectors for hitmaps
0138   // and hottest cell info
0139   template <typename T>
0140   double eHCALmatrix(const CaloGeometry* geo,
0141                      const HcalTopology* topology,
0142                      const DetId& det0,
0143                      edm::Handle<T>& hits,
0144                      int ieta,
0145                      int iphi,
0146                      int& nRecHits,
0147                      std::vector<int>& RH_ieta,
0148                      std::vector<int>& RH_iphi,
0149                      std::vector<double>& RH_ene,
0150                      GlobalPoint& gPosHotCell,
0151                      int useRaw) {
0152     bool debug = false;
0153     bool includeHO = false;
0154     HcalDetId hcid0(det0.rawId());
0155     HcalDetId hcid(hcid0.subdet(), hcid0.ieta(), hcid0.iphi(), 1);
0156     DetId det(hcid.rawId());
0157 
0158     std::vector<typename T::const_iterator> hit;
0159     spr::hitHCALmatrix(topology, det, hits, ieta, iphi, hit, includeHO, debug);
0160     double energyMax = -99.;
0161     int hottestIndex = -99;
0162 
0163     nRecHits = -99;
0164     double energySum = 0.0;
0165     std::set<int> uniqueIdset;
0166 
0167     for (unsigned int ihit = 0; ihit < hit.size(); ihit++) {
0168       // Sum energy getting hottest cell for later
0169       double energy = getEnergy(hit.at(ihit), useRaw);
0170       if (energy > energyMax) {
0171         energyMax = energy;
0172         hottestIndex = ihit;
0173       }
0174       energySum += energy;
0175       int eta(-99);
0176       int phi(-99);
0177       // Get eta phi for counting unique cells
0178       spr::getEtaPhi(hit[ihit], eta, phi);
0179       // Put eta phi in vectors
0180       RH_ieta.push_back(hit[ihit]->id().ieta());
0181       RH_iphi.push_back(hit[ihit]->id().iphi());
0182       RH_ene.push_back(hit[ihit]->energy());
0183       int uniqueEtaPhiId = 100 * eta + phi;
0184       uniqueIdset.insert(uniqueEtaPhiId);
0185     }
0186 
0187     // Get dist from center of cluster to hottest cell:
0188     if (hottestIndex != -99) {
0189       gPosHotCell = spr::getGpos(geo, hit.at(hottestIndex));
0190     }
0191 
0192     nRecHits = uniqueIdset.size();
0193     return energySum;
0194   }
0195 
0196   template <typename T>
0197   double eHCALmatrix(const CaloGeometry* geo,
0198                      const HcalTopology* topology,
0199                      const DetId& det0,
0200                      edm::Handle<T>& hits,
0201                      int ieta,
0202                      int iphi,
0203                      HcalDetId& hotCell,
0204                      bool includeHO,
0205                      int useRaw,
0206                      bool debug) {
0207     HcalDetId hcid0(det0.rawId());
0208     HcalDetId hcid(hcid0.subdet(), hcid0.ieta(), hcid0.iphi(), 1);
0209     DetId det(hcid.rawId());
0210 
0211     std::vector<typename T::const_iterator> hit;
0212     spr::hitHCALmatrix(topology, det, hits, ieta, iphi, hit, includeHO, debug);
0213 
0214     double energySum(0);
0215     for (unsigned int ihit = 0; ihit < hit.size(); ihit++)
0216       energySum += getEnergy(hit.at(ihit), useRaw);
0217 
0218     // Get hotCell ID
0219     hotCell = getHotCell(hit, includeHO, useRaw, debug);
0220     return energySum;
0221   }
0222 
0223   template <typename T>
0224   double energyHCALmatrixNew(const HcalTopology* topology,
0225                              const DetId& det,
0226                              edm::Handle<T>& hits,
0227                              int ieta,
0228                              int iphi,
0229                              bool includeHO,
0230                              double hbThr,
0231                              double heThr,
0232                              double hfThr,
0233                              double hoThr,
0234                              double tMin,
0235                              double tMax,
0236                              int useRaw,
0237                              bool debug) {
0238     std::vector<DetId> dets(1, det);
0239     std::vector<DetId> vdets = spr::matrixHCALIds(dets, topology, ieta, iphi, includeHO, false);
0240     if (debug) {
0241       edm::LogVerbatim("IsoTrack") << "matrixHCALIds::Total number of cells found is " << vdets.size();
0242       spr::debugHcalDets(0, vdets);
0243     }
0244     return spr::energyHCAL(vdets, hits, hbThr, heThr, hfThr, hoThr, tMin, tMax, useRaw, debug);
0245   }
0246 
0247   template <typename T>
0248   double energyHCALmatrixTotal(const HcalTopology* topology,
0249                                const DetId& det,
0250                                edm::Handle<T>& hits,
0251                                int ietaE,
0252                                int ietaW,
0253                                int iphiN,
0254                                int iphiS,
0255                                bool includeHO,
0256                                double hbThr,
0257                                double heThr,
0258                                double hfThr,
0259                                double hoThr,
0260                                double tMin,
0261                                double tMax,
0262                                int useRaw,
0263                                bool debug) {
0264     std::vector<DetId> dets(1, det);
0265     std::vector<DetId> vdets = spr::matrixHCALIds(dets, topology, ietaE, ietaW, iphiN, iphiS, includeHO, debug);
0266     return spr::energyHCAL(vdets, hits, hbThr, heThr, hfThr, hoThr, tMin, tMax, useRaw, debug);
0267   }
0268 
0269   template <typename T>
0270   void hitHCALmatrix(const HcalTopology* topology,
0271                      const DetId& det,
0272                      edm::Handle<T>& hits,
0273                      int ieta,
0274                      int iphi,
0275                      std::vector<typename T::const_iterator>& hitlist,
0276                      bool includeHO,
0277                      bool debug) {
0278     std::vector<DetId> dets(1, det);
0279     std::vector<DetId> vdets = spr::matrixHCALIds(dets, topology, ieta, iphi, includeHO, debug);
0280     spr::hitsHCAL(vdets, hits, hitlist, debug);
0281   }
0282 
0283   template <typename T>
0284   void hitHCALmatrixTotal(const HcalTopology* topology,
0285                           const DetId& det,
0286                           edm::Handle<T>& hits,
0287                           int ietaE,
0288                           int ietaW,
0289                           int iphiN,
0290                           int iphiS,
0291                           std::vector<typename T::const_iterator>& hitlist,
0292                           bool includeHO,
0293                           bool debug) {
0294     std::vector<DetId> dets(1, det);
0295     std::vector<DetId> vdets = spr::matrixHCALIds(dets, topology, ietaE, ietaW, iphiN, iphiS, includeHO, debug);
0296     spr::hitsHCAL(vdets, hits, hitlist, debug);
0297   }
0298 
0299   template <typename T>
0300   double energyHCAL(std::vector<DetId>& vdets,
0301                     edm::Handle<T>& hits,
0302                     double hbThr,
0303                     double heThr,
0304                     double hfThr,
0305                     double hoThr,
0306                     double tMin,
0307                     double tMax,
0308                     int useRaw,
0309                     bool debug) {
0310     int khit = 0;
0311     double energySum = 0;
0312     unsigned int maxsize = vdets.size();
0313     for (unsigned int i = 0; i < maxsize; i++) {
0314       std::vector<typename T::const_iterator> hit = spr::findHit(hits, vdets[i]);
0315       double energy = 0;
0316       int subdet = ((HcalDetId)(vdets[i].rawId())).subdet();
0317       double eThr = spr::eHCALThreshold(subdet, hbThr, heThr, hfThr, hoThr);
0318       for (unsigned int ihit = 0; ihit < hit.size(); ihit++) {
0319         if (hit[ihit] != hits->end()) {
0320           khit++;
0321           if (debug)
0322             edm::LogVerbatim("IsoTrack") << "energyHCAL:: Hit " << khit << " " << (HcalDetId)vdets[i] << " E "
0323                                          << hit[ihit]->energy() << " t " << hit[ihit]->time();
0324           if (hit[ihit]->time() > tMin && hit[ihit]->time() < tMax) {
0325             energy += getRawEnergy(hit[ihit], useRaw);
0326           }
0327         }
0328       }
0329       if (energy > eThr)
0330         energySum += energy;
0331     }
0332 
0333     if (debug)
0334       edm::LogVerbatim("IsoTrack") << "energyHCAL::Total Energy " << energySum << " from " << khit << " hits";
0335     return energySum;
0336   }
0337 
0338   template <typename T>
0339   void energyHCALCell(HcalDetId detID,
0340                       edm::Handle<T>& hits,
0341                       std::vector<std::pair<double, int> >& energyCell,
0342                       int maxDepth,
0343                       double hbThr,
0344                       double heThr,
0345                       double hfThr,
0346                       double hoThr,
0347                       double tMin,
0348                       double tMax,
0349                       int useRaw,
0350                       int depthHE,
0351                       bool debug) {
0352     energyCell.clear();
0353     int subdet = detID.subdet();
0354     double eThr = spr::eHCALThreshold(subdet, hbThr, heThr, hfThr, hoThr);
0355     bool hbhe = (detID.ietaAbs() == 16);
0356     if (debug)
0357       edm::LogVerbatim("IsoTrack") << "energyHCALCell: input ID " << detID << " MaxDepth " << maxDepth
0358                                    << " Threshold (E) " << eThr << " (T) " << tMin << ":" << tMax;
0359     for (int i = 0; i < maxDepth; i++) {
0360       HcalSubdetector subdet0 = (hbhe) ? ((i + 1 >= depthHE) ? HcalEndcap : HcalBarrel) : detID.subdet();
0361       HcalDetId hcid(subdet0, detID.ieta(), detID.iphi(), i + 1);
0362       DetId det(hcid.rawId());
0363       std::vector<typename T::const_iterator> hit = spr::findHit(hits, det);
0364       double energy(0);
0365       for (unsigned int ihit = 0; ihit < hit.size(); ++ihit) {
0366         if (hit[ihit]->time() > tMin && hit[ihit]->time() < tMax)
0367           energy += getRawEnergy(hit[ihit], useRaw);
0368         if (debug)
0369           edm::LogVerbatim("IsoTrack") << "energyHCALCell:: Hit[" << ihit << "] " << hcid << " E "
0370                                        << hit[ihit]->energy() << " t " << hit[ihit]->time();
0371       }
0372       if (debug)
0373         edm::LogVerbatim("IsoTrack") << "energyHCALCell:: Cell " << hcid << " E " << energy << " from " << hit.size()
0374                                      << " threshold " << eThr;
0375       if (energy > eThr && !(hit.empty())) {
0376         energyCell.push_back(std::pair<double, int>(energy, i + 1));
0377       }
0378     }
0379     if (debug) {
0380       edm::LogVerbatim("IsoTrack") << "energyHCALCell:: " << energyCell.size() << " entries from " << maxDepth
0381                                    << " depths:";
0382       std::ostringstream st1;
0383       for (unsigned int i = 0; i < energyCell.size(); ++i) {
0384         st1 << " [" << i << "] (" << energyCell[i].first << ":" << energyCell[i].second << ")";
0385       }
0386       edm::LogVerbatim("IsoTrack") << st1.str();
0387     }
0388   }
0389 
0390   template <typename T>
0391   void hitsHCAL(std::vector<DetId>& vdets,
0392                 edm::Handle<T>& hits,
0393                 std::vector<typename T::const_iterator>& hitlist,
0394                 bool debug) {
0395     unsigned int maxsize = vdets.size();
0396     for (unsigned int i = 0; i < maxsize; i++) {
0397       std::vector<typename T::const_iterator> hit = spr::findHit(hits, vdets[i]);
0398       hitlist.insert(hitlist.end(), hit.begin(), hit.end());
0399     }
0400 
0401     if (debug)
0402       edm::LogVerbatim("IsoTrack") << "hitsHCAL::Number of hits " << hitlist.size();
0403   }
0404 }  // namespace spr
[ Source navigation ] [ Diff markup ] [ Identifier search ] [ general search ]

This page was automatically generated by the 2.2.1 LXR engine.

The LXR team
Valid CSS 2.1! Valid HTML 4.01!