TrackIsoCalculator.cc

Go to the documentation of this file.

// ROOT includes
#include <Math/VectorUtil.h>

#include "RecoHI/HiEgammaAlgos/interface/TrackIsoCalculator.h"
#include "FWCore/MessageLogger/interface/MessageLogger.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"
#include "DataFormats/Math/interface/Vector3D.h"

using namespace edm;
using namespace reco;
using namespace std;

TrackIsoCalculator::TrackIsoCalculator(const edm::Event &iEvent,
                                       const edm::EventSetup &iSetup,
                                       const edm::Handle<reco::TrackCollection> trackLabel,
                                       const std::string trackQuality) {
  recCollection = trackLabel;
  trackQuality_ = trackQuality;
}

double TrackIsoCalculator::getTrackIso(const reco::Photon cluster,
                                       const double x,
                                       const double threshold,
                                       const double innerDR) {
  double TotalPt = 0;

  for (reco::TrackCollection::const_iterator recTrack = recCollection->begin(); recTrack != recCollection->end();
       recTrack++) {
    bool goodtrack = recTrack->quality(reco::TrackBase::qualityByName(trackQuality_));
    if (!goodtrack)
      continue;

    double pt = recTrack->pt();
    double dR2 = reco::deltaR2(cluster, *recTrack);
    if (dR2 >= (0.01 * x * x))
      continue;
    if (dR2 < innerDR * innerDR)
      continue;
    if (pt > threshold)
      TotalPt = TotalPt + pt;
  }

  return TotalPt;
}

double TrackIsoCalculator::getBkgSubTrackIso(const reco::Photon cluster,
                                             const double x,
                                             const double threshold,
                                             const double innerDR) {
  double SClusterEta = cluster.eta();
  double TotalPt = 0;

  TotalPt = 0;

  for (reco::TrackCollection::const_iterator recTrack = recCollection->begin(); recTrack != recCollection->end();
       recTrack++) {
    bool goodtrack = recTrack->quality(reco::TrackBase::qualityByName(trackQuality_));
    if (!goodtrack)
      continue;

    double pt = recTrack->pt();
    double eta2 = recTrack->eta();
    double dEta = fabs(eta2 - SClusterEta);
    double dR2 = reco::deltaR2(cluster, *recTrack);
    if (dEta >= 0.1 * x)
      continue;
    if (dR2 < innerDR * innerDR)
      continue;

    if (pt > threshold)
      TotalPt = TotalPt + pt;
  }

  double Tx = getTrackIso(cluster, x, threshold, innerDR);
  double CTx = (Tx - TotalPt / 40.0 * x) * (1 / (1 - x / 40.));

  return CTx;
}