EcalPFClusterIsolation.cc

Go to the documentation of this file.

//*****************************************************************************
// File:      EgammaRecHitIsolation.cc
// ----------------------------------------------------------------------------
// OrigAuth:  Matthias Mozer, hacked by Sam Harper (ie the ugly stuff is mine)
// Institute: IIHE-VUB, RAL
//=============================================================================
//*****************************************************************************

#include "RecoEgamma/EgammaIsolationAlgos/interface/EcalPFClusterIsolation.h"

#include "DataFormats/RecoCandidate/interface/RecoEcalCandidate.h"
#include "DataFormats/RecoCandidate/interface/RecoChargedCandidate.h"
#include "DataFormats/EgammaCandidates/interface/GsfElectron.h"
#include "DataFormats/EgammaCandidates/interface/Photon.h"

#include <DataFormats/Math/interface/deltaR.h>

template <typename T1>
EcalPFClusterIsolation<T1>::EcalPFClusterIsolation(double drMax,
                                                   double drVetoBarrel,
                                                   double drVetoEndcap,
                                                   double etaStripBarrel,
                                                   double etaStripEndcap,
                                                   double energyBarrel,
                                                   double energyEndcap)
    : drMax_(drMax),
      drVetoBarrel_(drVetoBarrel),
      drVetoEndcap_(drVetoEndcap),
      etaStripBarrel_(etaStripBarrel),
      etaStripEndcap_(etaStripEndcap),
      energyBarrel_(energyBarrel),
      energyEndcap_(energyEndcap) {}

template <typename T1>
EcalPFClusterIsolation<T1>::~EcalPFClusterIsolation() {}

template <typename T1>
double EcalPFClusterIsolation<T1>::getSum(const T1Ref candRef, edm::Handle<reco::PFClusterCollection> clusterHandle) {
  drVeto2_ = -1.;
  float etaStrip = -1;

  if (fabs(candRef->eta()) < 1.479) {
    drVeto2_ = drVetoBarrel_ * drVetoBarrel_;
    etaStrip = etaStripBarrel_;
  } else {
    drVeto2_ = drVetoEndcap_ * drVetoEndcap_;
    etaStrip = etaStripEndcap_;
  }

  float etSum = 0;
  for (size_t i = 0; i < clusterHandle->size(); i++) {
    reco::PFClusterRef pfclu(clusterHandle, i);

    if (fabs(candRef->eta()) < 1.479) {
      if (fabs(pfclu->pt()) < energyBarrel_)
        continue;
    } else {
      if (fabs(pfclu->energy()) < energyEndcap_)
        continue;
    }

    float dEta = fabs(candRef->eta() - pfclu->eta());
    if (dEta < etaStrip)
      continue;
    if (not computedRVeto(candRef, pfclu))
      continue;

    etSum += pfclu->pt();
  }

  return etSum;
}

template <typename T1>
bool EcalPFClusterIsolation<T1>::computedRVeto(T1Ref candRef, reco::PFClusterRef pfclu) {
  float dR2 = deltaR2(candRef->eta(), candRef->phi(), pfclu->eta(), pfclu->phi());
  if (dR2 > (drMax_ * drMax_))
    return false;

  if (candRef->superCluster().isNonnull()) {
    // Exclude clusters that are part of the candidate
    for (reco::CaloCluster_iterator it = candRef->superCluster()->clustersBegin();
         it != candRef->superCluster()->clustersEnd();
         ++it) {
      if ((*it)->seed() == pfclu->seed()) {
        return false;
      }
    }
  }

  return true;
}

template <>
bool EcalPFClusterIsolation<reco::RecoChargedCandidate>::computedRVeto(T1Ref candRef, reco::PFClusterRef pfclu) {
  float dR2 = deltaR2(candRef->eta(), candRef->phi(), pfclu->eta(), pfclu->phi());
  if (dR2 > (drMax_ * drMax_) || dR2 < drVeto2_)
    return false;
  else
    return true;
}

template class EcalPFClusterIsolation<reco::RecoEcalCandidate>;
template class EcalPFClusterIsolation<reco::RecoChargedCandidate>;
template class EcalPFClusterIsolation<reco::Photon>;
template class EcalPFClusterIsolation<reco::GsfElectron>;