HiEgammaAlgos/src/HcalRechitIsoCalculator.cc

0001 #include "RecoHI/HiEgammaAlgos/interface/HcalRechitIsoCalculator.h"
0002
0003 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0004
0005 #include "Geometry/CaloGeometry/interface/CaloSubdetectorGeometry.h"
0006
0007 #include "DataFormats/Common/interface/Handle.h"
0008 #include "DataFormats/EgammaReco/interface/BasicCluster.h"
0009 #include "DataFormats/EgammaReco/interface/SuperCluster.h"
0010
0011 using namespace edm;
0012 using namespace reco;
0013
0014 HcalRechitIsoCalculator::HcalRechitIsoCalculator(const CaloGeometry *geometry,
0015                                                  const edm::Handle<HBHERecHitCollection> hbhe,
0016                                                  const edm::Handle<HFRecHitCollection> hf,
0017                                                  const edm::Handle<HORecHitCollection> ho) {
0018   if (hf.isValid())
0019     fHFRecHits_ = hf.product();
0020   else
0021     fHFRecHits_ = nullptr;
0022
0023   if (ho.isValid())
0024     fHORecHits_ = ho.product();
0025   else
0026     fHORecHits_ = nullptr;
0027
0028   if (hbhe.isValid())
0029     fHBHERecHits_ = hbhe.product();
0030   else
0031     fHBHERecHits_ = nullptr;
0032
0033   geometry_ = geometry;
0034 }
0035
0036 double HcalRechitIsoCalculator::getHcalRechitIso(const reco::SuperClusterRef cluster,
0037                                                  const double x,
0038                                                  const double threshold,
0039                                                  const double innerR) {
0040   if (!fHBHERecHits_) {
0041     return -100;
0042   }
0043
0044   double TotalEt = 0;
0045
0046   for (size_t index = 0; index < fHBHERecHits_->size(); index++) {
0047     const HBHERecHit &rechit = (*fHBHERecHits_)[index];
0048     const DetId &detid = rechit.id();
0049     const GlobalPoint &hitpoint = geometry_->getPosition(detid);
0050     double eta = hitpoint.eta();
0051
0052     double dR2 = reco::deltaR2(*cluster, hitpoint);
0053     // veto inner cone///////////////
0054     if (dR2 < innerR * innerR)
0055       continue;
0056     /////////////////////////////////
0057     if (dR2 < (x * x * 0.01)) {
0058       double et = rechit.energy() / cosh(eta);
0059       if (et < threshold)
0060         et = 0;
0061       TotalEt += et;
0062     }
0063   }
0064
0065   return TotalEt;
0066 }
0067
0068 double HcalRechitIsoCalculator::getBkgSubHcalRechitIso(const reco::SuperClusterRef cluster,
0069                                                        const double x,
0070                                                        const double threshold,
0071                                                        const double innerR) {
0072   if (!fHBHERecHits_) {
0073     return -100;
0074   }
0075
0076   double SClusterEta = cluster->eta();
0077   double TotalEt = 0;
0078
0079   for (size_t index = 0; index < fHBHERecHits_->size(); index++) {
0080     const HBHERecHit &rechit = (*fHBHERecHits_)[index];
0081     const DetId &detid = rechit.id();
0082     const GlobalPoint &hitpoint = geometry_->getPosition(detid);
0083     double eta = hitpoint.eta();
0084     double dEta = fabs(eta - SClusterEta);
0085
0086     if (dEta < x * 0.1) {
0087       double et = rechit.energy() / cosh(eta);
0088       if (et < threshold)
0089         et = 0;
0090       TotalEt += et;
0091     }
0092   }
0093
0094   double Rx = getHcalRechitIso(cluster, x, threshold, innerR);
0095   double CRx = (Rx - TotalEt * (0.01 * x * x - innerR * innerR) / (2 * 2 * 0.1 * x)) * (1 / (1 - x / 40.));
0096
0097   return CRx;
0098 }