HcalRechitIsoCalculator.cc

Go to the documentation of this file.

#include "RecoHI/HiEgammaAlgos/interface/HcalRechitIsoCalculator.h"

#include "FWCore/MessageLogger/interface/MessageLogger.h"

#include "Geometry/Records/interface/IdealGeometryRecord.h"
#include "Geometry/Records/interface/CaloGeometryRecord.h"
#include "Geometry/CaloGeometry/interface/CaloSubdetectorGeometry.h"

#include "DataFormats/Common/interface/Handle.h"
#include "DataFormats/EgammaReco/interface/BasicCluster.h"
#include "DataFormats/EgammaReco/interface/SuperCluster.h"


using namespace edm;
using namespace reco;

HcalRechitIsoCalculator::HcalRechitIsoCalculator(const edm::Event &iEvent, const edm::EventSetup &iSetup, const edm::Handle<HBHERecHitCollection> hbhe, const edm::Handle<HFRecHitCollection> hf, const edm::Handle<HORecHitCollection> ho)
{
  if(hf.isValid())
    fHFRecHits_ = hf.product();
  else
    fHFRecHits_ = nullptr;

  if(ho.isValid())
    fHORecHits_ = ho.product();
  else
    fHORecHits_ = nullptr;

  if(hbhe.isValid())
    fHBHERecHits_ = hbhe.product();
  else
    fHBHERecHits_ = nullptr;

  ESHandle<CaloGeometry> geometryHandle;
  iSetup.get<CaloGeometryRecord>().get(geometryHandle);
  if(geometryHandle.isValid())
    geometry_ = geometryHandle.product();
  else
    geometry_ = nullptr;

}


double HcalRechitIsoCalculator::getHcalRechitIso(const reco::SuperClusterRef cluster, const double x, const double threshold, const double innerR )
{
   if(!fHBHERecHits_) {
     return -100;
   }

   double TotalEt = 0;

   for(size_t index = 0; index < fHBHERecHits_->size(); index++) {
      const HBHERecHit &rechit = (*fHBHERecHits_)[index];
      const DetId &detid = rechit.id();
      const GlobalPoint& hitpoint = geometry_->getPosition(detid);
      double eta = hitpoint.eta();

      double dR2 = reco::deltaR2(*cluster, hitpoint);
      // veto inner cone///////////////
      if ( dR2 < innerR*innerR )  continue;
      if (dR2< (x*x*0.01)) {
     double et = rechit.energy()/cosh(eta);
         if (et<threshold) et=0;
         TotalEt += et;
      }
   }

   return TotalEt;
}

double HcalRechitIsoCalculator::getBkgSubHcalRechitIso(const reco::SuperClusterRef cluster, const double x, const double threshold, const double innerR )
{
   if(!fHBHERecHits_) {
      return -100;
   }

   double SClusterEta = cluster->eta();
   double TotalEt = 0;

   for(size_t index = 0; index < fHBHERecHits_->size(); index++) {
      const HBHERecHit &rechit = (*fHBHERecHits_)[index];
      const DetId &detid = rechit.id();
      const GlobalPoint& hitpoint = geometry_->getPosition(detid);
      double eta = hitpoint.eta();
      double dEta = fabs(eta-SClusterEta);

      if (dEta<x*0.1) {
         double et = rechit.energy()/cosh(eta);
         if (et<threshold) et=0;
         TotalEt += et;
      }
   }

   double Rx = getHcalRechitIso(cluster,x,threshold,innerR);
   double CRx = (Rx - TotalEt * (0.01*x*x - innerR*innerR) / (2 * 2 * 0.1 * x))*(1/(1-x/40.)) ;

   return CRx;
}