Back to home page

Project CMSSW displayed by LXR
 
 

    

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

0001 #include "Calibration/IsolatedParticles/interface/FindCaloHit.h"
0002 #include "Calibration/IsolatedParticles/interface/FindCaloHitCone.h"
0003 #include "Calibration/IsolatedParticles/interface/FindDistCone.h"
0004 #include "DataFormats/HcalDetId/interface/HcalSubdetector.h"
0005 #include "DataFormats/HcalDetId/interface/HcalDetId.h"
0006 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0007 #include <iostream>
0008 
0009 namespace spr {
0010 
0011   // Basic cone energy cluster for hcal simhits and hcal rechits
0012   template <typename T>
0013   double eCone_hcal(const CaloGeometry* geo,
0014                     edm::Handle<T>& hits,
0015                     const GlobalPoint& hpoint1,
0016                     const GlobalPoint& point1,
0017                     double dR,
0018                     const GlobalVector& trackMom,
0019                     int& nRecHits,
0020                     double hbThr,
0021                     double heThr,
0022                     double hfThr,
0023                     double hoThr,
0024                     double tMin,
0025                     double tMax,
0026                     int detOnly,
0027                     int useRaw,
0028                     bool debug) {
0029     std::vector<typename T::const_iterator> hit = spr::findHitCone(geo, hits, hpoint1, point1, dR, trackMom, debug);
0030 
0031     nRecHits = -99;
0032     double energySum = 0.0;
0033     std::set<int> uniqueIdset;
0034     for (auto const& ihit : hit) {
0035       int subdet = DetId(ihit->id()).subdetId();
0036       if (detOnly < 0 || subdet == detOnly) {
0037         double eTower = spr::getEnergy(ihit, useRaw);
0038         double tt = ihit->time();
0039         bool ok = (eTower > hbThr);
0040         if (subdet == (int)(HcalEndcap))
0041           ok = (eTower > heThr);
0042         else if (subdet == (int)(HcalForward))
0043           ok = (eTower > hfThr);
0044         else if (subdet == (int)(HcalOuter))
0045           ok = (eTower > hoThr);
0046         if (ok && tt > tMin && tt < tMax) {
0047           energySum += eTower;
0048           int eta(-99);
0049           int phi(-99);
0050           spr::getEtaPhi(ihit, eta, phi);
0051           int uniqueEtaPhiId = 100 * eta + phi;
0052           uniqueIdset.insert(uniqueEtaPhiId);
0053         }
0054       }
0055     }
0056     nRecHits = uniqueIdset.size();
0057     if (debug) {
0058       edm::LogVerbatim("IsoTrack") << "eCone_hcal: Energy " << energySum << " from " << nRecHits << " in subdet "
0059                                    << detOnly;
0060     }
0061     return energySum;
0062   }
0063 
0064   // Cone energy cluster for hcal simhits and hcal rechits
0065   // that returns vector of rechit IDs and hottest cell info
0066   template <typename T>
0067   double eCone_hcal(const CaloGeometry* geo,
0068                     edm::Handle<T>& hits,
0069                     const GlobalPoint& hpoint1,
0070                     const GlobalPoint& point1,
0071                     double dR,
0072                     const GlobalVector& trackMom,
0073                     int& nRecHits,
0074                     std::vector<DetId>& coneRecHitDetIds,
0075                     double& distFromHotCell,
0076                     int& ietaHotCell,
0077                     int& iphiHotCell,
0078                     GlobalPoint& gposHotCell,
0079                     int detOnly,
0080                     int useRaw,
0081                     bool debug) {
0082     std::vector<typename T::const_iterator> hit = spr::findHitCone(geo, hits, hpoint1, point1, dR, trackMom, debug);
0083 
0084     double energyMax = -99.;
0085     int hottestIndex = -99;
0086     ietaHotCell = -99;
0087     iphiHotCell = -99;
0088     distFromHotCell = -99.;
0089 
0090     nRecHits = -99;
0091     double energySum = 0.0;
0092     std::set<int> uniqueIdset;
0093     for (unsigned int ihit = 0; ihit < hit.size(); ihit++) {
0094       int subdet = DetId(hit[ihit]->id()).subdetId();
0095       if (detOnly < 0 || subdet == detOnly) {
0096         // Sum energy getting hottest cell for later
0097         double energy = spr::getEnergy(hit.at(ihit), useRaw);
0098         if (energy > energyMax) {
0099           energyMax = energy;
0100           hottestIndex = ihit;
0101         }
0102         energySum += energy;
0103         // Don't double count depth in rechit multiplicity
0104         int eta(-99);
0105         int phi(-99);
0106         spr::getEtaPhi(hit[ihit], eta, phi);
0107         int uniqueEtaPhiId = 100 * eta + phi;
0108         uniqueIdset.insert(uniqueEtaPhiId);
0109         // Get list of det ids
0110         coneRecHitDetIds.emplace_back(hit[ihit]->id());
0111       }
0112     }
0113 
0114     // Get dist from center of cluster to hottest cell:
0115     if (hottestIndex != -99) {
0116       gposHotCell = spr::getGpos(geo, hit.at(hottestIndex));
0117       ietaHotCell = hit.at(hottestIndex)->id().ieta();
0118       iphiHotCell = hit.at(hottestIndex)->id().iphi();
0119       distFromHotCell = spr::getDistInPlaneTrackDir(hpoint1, trackMom, gposHotCell);
0120     }
0121 
0122     nRecHits = uniqueIdset.size();
0123     if (debug) {
0124       edm::LogVerbatim("IsoTrack") << "eCone_hcal: Energy " << energySum << " from " << nRecHits << ":"
0125                                    << coneRecHitDetIds.size() << " in subdet " << detOnly;
0126       edm::LogVerbatim("IsoTrack") << "HotCell " << ietaHotCell << ":" << iphiHotCell << " dist " << distFromHotCell;
0127     }
0128     return energySum;
0129   }
0130 
0131   template <typename T>
0132   double eCone_hcal(const CaloGeometry* geo,
0133                     edm::Handle<T>& hits,
0134                     const GlobalPoint& hpoint1,
0135                     const GlobalPoint& point1,
0136                     double dR,
0137                     const GlobalVector& trackMom,
0138                     int& nRecHits,
0139                     std::vector<DetId>& coneRecHitDetIds,
0140                     std::vector<double>& eHit,
0141                     int useRaw,
0142                     bool debug) {
0143     std::vector<typename T::const_iterator> hit = spr::findHitCone(geo, hits, hpoint1, point1, dR, trackMom, debug);
0144     nRecHits = -99;
0145     double energySum = 0.0;
0146     std::set<int> uniqueIdset;
0147     for (auto const& ihit : hit) {
0148       // Sum energy getting hottest cell for later
0149       double energy = spr::getEnergy(ihit, useRaw);
0150       energySum += energy;
0151       // Don't double count depth in rechit multiplicity
0152       int eta(-99), phi(-99);
0153       spr::getEtaPhi(ihit, eta, phi);
0154       int uniqueEtaPhiId = 100 * eta + phi;
0155       uniqueIdset.insert(uniqueEtaPhiId);
0156       // Get list of det ids
0157       coneRecHitDetIds.emplace_back(ihit->id());
0158       eHit.emplace_back(energy);
0159     }
0160 
0161     nRecHits = uniqueIdset.size();
0162     if (debug) {
0163       edm::LogVerbatim("IsoTrack") << "eCone_hcal: Energy " << energySum << " from " << nRecHits << ":"
0164                                    << coneRecHitDetIds.size();
0165       for (unsigned int k = 0; k < eHit.size(); ++k)
0166         edm::LogVerbatim("IsoTrack") << "Hit[" << k << "] " << HcalDetId(coneRecHitDetIds[k]) << " energy " << eHit[k];
0167     }
0168     return energySum;
0169   }
0170 
0171   // Cone energy cluster for hcal simhits and hcal rechits
0172   // that returns vector of rechit IDs and hottest cell info
0173   // AND info for making "hit maps"
0174   template <typename T>
0175   double eCone_hcal(const CaloGeometry* geo,
0176                     edm::Handle<T>& hits,
0177                     const GlobalPoint& hpoint1,
0178                     const GlobalPoint& point1,
0179                     double dR,
0180                     const GlobalVector& trackMom,
0181                     int& nRecHits,
0182                     std::vector<int>& RH_ieta,
0183                     std::vector<int>& RH_iphi,
0184                     std::vector<double>& RH_ene,
0185                     std::vector<DetId>& coneRecHitDetIds,
0186                     double& distFromHotCell,
0187                     int& ietaHotCell,
0188                     int& iphiHotCell,
0189                     GlobalPoint& gposHotCell,
0190                     int detOnly,
0191                     int useRaw,
0192                     bool debug) {
0193     std::vector<typename T::const_iterator> hit = spr::findHitCone(geo, hits, hpoint1, point1, dR, trackMom, debug);
0194 
0195     double energyMax = -99.;
0196     int hottestIndex = -99;
0197     ietaHotCell = -99;
0198     iphiHotCell = -99;
0199     distFromHotCell = -99.;
0200 
0201     nRecHits = -99;
0202     double energySum = 0.0;
0203     std::set<int> uniqueIdset;
0204     for (unsigned int ihit = 0; ihit < hit.size(); ihit++) {
0205       int subdet = DetId(hit[ihit]->id()).subdetId();
0206       if (detOnly < 0 || subdet == detOnly) {
0207         // Sum energy getting hottest cell for later
0208         double energy = spr::getEnergy(hit.at(ihit), useRaw);
0209         if (energy > energyMax) {
0210           energyMax = energy;
0211           hottestIndex = ihit;
0212         }
0213         energySum += energy;
0214 
0215         // Get eta phi for counting unique cells
0216         int eta(-99);
0217         int phi(-99);
0218         spr::getEtaPhi(hit[ihit], eta, phi);
0219         // Put eta phi in vectors
0220         spr::getEtaPhi(hit[ihit], RH_ieta, RH_iphi, RH_ene);
0221         int uniqueEtaPhiId = 100 * eta + phi;
0222         uniqueIdset.insert(uniqueEtaPhiId);
0223         // Get vector of detids.
0224         coneRecHitDetIds.emplace_back(hit[ihit]->id());
0225       }
0226     }
0227     nRecHits = uniqueIdset.size();
0228     // Get dist from center of cluster to hottest cell:
0229     if (hottestIndex != -99) {
0230       gposHotCell = spr::getGpos(geo, hit.at(hottestIndex));
0231       ietaHotCell = hit.at(hottestIndex)->id().ieta();
0232       iphiHotCell = hit.at(hottestIndex)->id().iphi();
0233       distFromHotCell = spr::getDistInPlaneTrackDir(hpoint1, trackMom, gposHotCell);
0234     }
0235 
0236     nRecHits = uniqueIdset.size();
0237     if (debug) {
0238       edm::LogVerbatim("IsoTrack") << "eCone_hcal: Energy " << energySum << " from " << nRecHits << ":"
0239                                    << coneRecHitDetIds.size() << " in subdet " << detOnly;
0240       edm::LogVerbatim("IsoTrack") << "HotCell " << ietaHotCell << ":" << iphiHotCell << " dist " << distFromHotCell;
0241     }
0242     return energySum;
0243   }
0244 
0245   //
0246   template <typename T>
0247   double eCone_ecal(const CaloGeometry* geo,
0248                     edm::Handle<T>& barrelhits,
0249                     edm::Handle<T>& endcaphits,
0250                     const GlobalPoint& hpoint1,
0251                     const GlobalPoint& point1,
0252                     double dR,
0253                     const GlobalVector& trackMom,
0254                     int& nRecHits,
0255                     double ebThr,
0256                     double eeThr,
0257                     double tMin,
0258                     double tMax,
0259                     bool debug) {
0260     std::vector<typename T::const_iterator> hit =
0261         spr::findHitCone(geo, barrelhits, endcaphits, hpoint1, point1, dR, trackMom, debug);
0262 
0263     // No depth in Ecal so nRecHits = hit.size()
0264     nRecHits = hit.size();
0265     double energySum = 0.0;
0266     for (auto const& ihit : hit) {
0267       bool ok = true;
0268       double eTower = ihit->energy();
0269       double tt = ihit->time();
0270       if (DetId(ihit->id()).subdetId() == EcalBarrel)
0271         ok = (eTower > ebThr);
0272       else if (DetId(ihit->id()).subdetId() == EcalEndcap)
0273         ok = (eTower > eeThr);
0274       // No need for hcal sampling weight
0275       if (ok && tt > tMin && tt < tMax)
0276         energySum += eTower;
0277     }
0278 
0279     if (debug) {
0280       edm::LogVerbatim("IsoTrack") << "eCone_ecal: Energy " << energySum << " from " << hit.size() << " hits";
0281     }
0282     return energySum;
0283   }
0284 
0285   //
0286   template <typename T>
0287   double eCone_ecal(const CaloGeometry* geo,
0288                     edm::Handle<T>& barrelhits,
0289                     edm::Handle<T>& endcaphits,
0290                     const GlobalPoint& hpoint1,
0291                     const GlobalPoint& point1,
0292                     double dR,
0293                     const GlobalVector& trackMom,
0294                     std::vector<DetId>& coneRecHitDetIds,
0295                     std::vector<double>& eHit,
0296                     double ebThr,
0297                     double eeThr,
0298                     double tMin,
0299                     double tMax,
0300                     bool debug) {
0301     std::vector<typename T::const_iterator> hit =
0302         spr::findHitCone(geo, barrelhits, endcaphits, hpoint1, point1, dR, trackMom, debug);
0303 
0304     // No depth in Ecal
0305     double energySum = 0.0;
0306     for (auto const& ihit : hit) {
0307       bool ok = true;
0308       double eTower = ihit->energy();
0309       double tt = ihit->time();
0310       if (DetId(ihit->id()).subdetId() == EcalBarrel)
0311         ok = (eTower > ebThr);
0312       else if (DetId(ihit->id()).subdetId() == EcalEndcap)
0313         ok = (eTower > eeThr);
0314       // No need for hcal sampling weight
0315       if (tt > tMin && tt < tMax) {
0316         if (ok)
0317           energySum += eTower;
0318         // Get list of det ids
0319         coneRecHitDetIds.emplace_back(ihit->id());
0320         eHit.emplace_back(eTower);
0321       }
0322     }
0323 
0324     if (debug) {
0325       edm::LogVerbatim("IsoTrack") << "eCone_ecal: Energy " << energySum << " from " << hit.size() << " hits";
0326     }
0327     return energySum;
0328   }
0329 
0330 }  // 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!