IsolatedParticles/src/FindCaloHitCone.cc

0001 #include "Calibration/IsolatedParticles/interface/CaloConstants.h"
0002 #include "Calibration/IsolatedParticles/interface/FindCaloHitCone.h"
0003 #include "Calibration/IsolatedParticles/interface/FindDistCone.h"
0004 #include "DataFormats/EcalDetId/interface/EcalSubdetector.h"
0005 #include "Geometry/HcalTowerAlgo/interface/HcalGeometry.h"
0006 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0007 #include <iostream>
0008
0009 namespace spr {
0010
0011   //Ecal Endcap OR Barrel RecHits
0012   std::vector<EcalRecHitCollection::const_iterator> findCone(const CaloGeometry* geo,
0013                                                              edm::Handle<EcalRecHitCollection>& hits,
0014                                                              const GlobalPoint& hpoint1,
0015                                                              const GlobalPoint& point1,
0016                                                              double dR,
0017                                                              const GlobalVector& trackMom,
0018                                                              bool debug) {
0019     std::vector<EcalRecHitCollection::const_iterator> hit;
0020
0021     for (EcalRecHitCollection::const_iterator j = hits->begin(); j != hits->end(); j++) {
0022       bool keepHit = false;
0023
0024       if (j->id().subdetId() == EcalEndcap) {
0025         EEDetId EEid = EEDetId(j->id());
0026         const GlobalPoint& rechitPoint = geo->getPosition(EEid);
0027         if (spr::getDistInPlaneTrackDir(point1, trackMom, rechitPoint, debug) < dR)
0028           keepHit = true;
0029       } else if (j->id().subdetId() == EcalBarrel) {
0030         EBDetId EBid = EBDetId(j->id());
0031         const GlobalPoint& rechitPoint = geo->getPosition(EBid);
0032         if (spr::getDistInPlaneTrackDir(point1, trackMom, rechitPoint, debug) < dR)
0033           keepHit = true;
0034       }
0035
0036       if (keepHit)
0037         hit.push_back(j);
0038     }
0039     return hit;
0040   }
0041
0042   // Ecal Endcap AND Barrel RecHits
0043   std::vector<EcalRecHitCollection::const_iterator> findCone(const CaloGeometry* geo,
0044                                                              edm::Handle<EcalRecHitCollection>& barrelhits,
0045                                                              edm::Handle<EcalRecHitCollection>& endcaphits,
0046                                                              const GlobalPoint& hpoint1,
0047                                                              const GlobalPoint& point1,
0048                                                              double dR,
0049                                                              const GlobalVector& trackMom,
0050                                                              bool debug) {
0051     std::vector<EcalRecHitCollection::const_iterator> hit;
0052
0053     // Only check both barrel and endcap when track is near transition
0054     // region: 1.479-2*0.087 < trkEta < 1.479+2*0.087
0055
0056     bool doBarrel = false, doEndcap = false;
0057     if (std::abs(point1.eta()) < (spr::etaBEEcal + 2 * spr::deltaEta))
0058       doBarrel = true;  // 1.479+2*0.087
0059     if (std::abs(point1.eta()) > (spr::etaBEEcal - 2 * spr::deltaEta))
0060       doEndcap = true;  // 1.479-2*0.087
0061
0062     if (doBarrel) {
0063       for (EcalRecHitCollection::const_iterator j = barrelhits->begin(); j != barrelhits->end(); j++) {
0064         bool keepHit = false;
0065         if (j->id().subdetId() == EcalBarrel) {
0066           EBDetId EBid = EBDetId(j->id());
0067           const GlobalPoint& rechitPoint = geo->getPosition(EBid);
0068           if (spr::getDistInPlaneTrackDir(point1, trackMom, rechitPoint, debug) < dR)
0069             keepHit = true;
0070         } else {
0071           edm::LogWarning("IsoTrack") << "PROBLEM : Endcap RecHits in Barrel Collection!?";
0072         }
0073         if (keepHit)
0074           hit.push_back(j);
0075       }
0076     }  // doBarrel
0077
0078     if (doEndcap) {
0079       for (EcalRecHitCollection::const_iterator j = endcaphits->begin(); j != endcaphits->end(); j++) {
0080         bool keepHit = false;
0081
0082         if (j->id().subdetId() == EcalEndcap) {
0083           EEDetId EEid = EEDetId(j->id());
0084           const GlobalPoint& rechitPoint = geo->getPosition(EEid);
0085           if (spr::getDistInPlaneTrackDir(point1, trackMom, rechitPoint, debug) < dR)
0086             keepHit = true;
0087         } else {
0088           edm::LogWarning("IsoTrack") << "PROBLEM : Barrel RecHits in Endcap Collection!?";
0089         }
0090         if (keepHit)
0091           hit.push_back(j);
0092       }
0093     }  // doEndcap
0094
0095     return hit;
0096   }
0097
0098   //HBHE RecHits
0099   std::vector<HBHERecHitCollection::const_iterator> findCone(const CaloGeometry* geo,
0100                                                              edm::Handle<HBHERecHitCollection>& hits,
0101                                                              const GlobalPoint& hpoint1,
0102                                                              const GlobalPoint& point1,
0103                                                              double dR,
0104                                                              const GlobalVector& trackMom,
0105                                                              bool debug) {
0106     std::vector<HBHERecHitCollection::const_iterator> hit;
0107     // Loop over Hcal RecHits
0108     for (HBHERecHitCollection::const_iterator j = hits->begin(); j != hits->end(); j++) {
0109       DetId detId(j->id());
0110       const GlobalPoint rechitPoint =
0111           (static_cast<const HcalGeometry*>(geo->getSubdetectorGeometry(detId)))->getPosition(detId);
0112       if (spr::getDistInPlaneTrackDir(hpoint1, trackMom, rechitPoint, debug) < dR)
0113         hit.push_back(j);
0114     }
0115     return hit;
0116   }
0117
0118   // PCalo SimHits
0119   std::vector<edm::PCaloHitContainer::const_iterator> findCone(const CaloGeometry* geo,
0120                                                                edm::Handle<edm::PCaloHitContainer>& hits,
0121                                                                const GlobalPoint& hpoint1,
0122                                                                const GlobalPoint& point1,
0123                                                                double dR,
0124                                                                const GlobalVector& trackMom,
0125                                                                bool debug) {
0126     std::vector<edm::PCaloHitContainer::const_iterator> hit;
0127     edm::PCaloHitContainer::const_iterator ihit;
0128     for (ihit = hits->begin(); ihit != hits->end(); ihit++) {
0129       DetId detId(ihit->id());
0130       const GlobalPoint rechitPoint =
0131           (detId.det() == DetId::Hcal)
0132               ? (static_cast<const HcalGeometry*>(geo->getSubdetectorGeometry(detId)))->getPosition(detId)
0133               : geo->getPosition(detId);
0134       if (spr::getDistInPlaneTrackDir(hpoint1, trackMom, rechitPoint, debug) < dR) {
0135         hit.push_back(ihit);
0136       }
0137     }
0138     return hit;
0139   }
0140
0141 }  // namespace spr