EgammaTrackExtractor.cc

Go to the documentation of this file.

#include "RecoEgamma/EgammaIsolationAlgos/plugins/EgammaTrackExtractor.h"

#include "RecoEgamma/EgammaIsolationAlgos/interface/EgammaRange.h"
#include "RecoEgamma/EgammaIsolationAlgos/interface/EgammaRecHitIsolation.h"
#include "DataFormats/RecoCandidate/interface/IsoDepositDirection.h"
#include "RecoEgamma/EgammaIsolationAlgos/interface/EgammaTrackSelector.h"
#include "DataFormats/Common/interface/Handle.h"
#include "DataFormats/TrackReco/interface/Track.h"
#include "DataFormats/TrackReco/interface/TrackFwd.h"
#include "FWCore/MessageLogger/interface/MessageLogger.h"
#include "DataFormats/BeamSpot/interface/BeamSpot.h"
#include "Geometry/Records/interface/CaloGeometryRecord.h"
#include "DataFormats/EgammaReco/interface/SuperCluster.h"
#include "DataFormats/EgammaReco/interface/SuperClusterFwd.h"
#include "DataFormats/EgammaCandidates/interface/GsfElectron.h"
#include "DataFormats/EgammaCandidates/interface/GsfElectronFwd.h"
#include "DataFormats/GsfTrackReco/interface/GsfTrack.h"
#include "DataFormats/GsfTrackReco/interface/GsfTrackFwd.h"
#include "DataFormats/RecoCandidate/interface/RecoCandidate.h"
#include "CommonTools/Statistics/interface/ChiSquaredProbability.h"

using namespace edm;
using namespace std;
using namespace reco;
using namespace egammaisolation;
using reco::isodeposit::Direction;

EgammaTrackExtractor::EgammaTrackExtractor(const ParameterSet& par, edm::ConsumesCollector& iC)
    : theTrackCollectionToken(iC.consumes<View<Track> >(par.getParameter<edm::InputTag>("inputTrackCollection"))),
      theDepositLabel(par.getUntrackedParameter<std::string>("DepositLabel")),
      theDiff_r(par.getParameter<double>("Diff_r")),
      theDiff_z(par.getParameter<double>("Diff_z")),
      theDR_Max(par.getParameter<double>("DR_Max")),
      theDR_Veto(par.getParameter<double>("DR_Veto")),
      theBeamlineOption(par.getParameter<std::string>("BeamlineOption")),
      theBeamSpotToken(iC.mayConsume<reco::BeamSpot>(par.getParameter<edm::InputTag>("BeamSpotLabel"))),
      theNHits_Min(par.getParameter<unsigned int>("NHits_Min")),
      theChi2Ndof_Max(par.getParameter<double>("Chi2Ndof_Max")),
      theChi2Prob_Min(par.getParameter<double>("Chi2Prob_Min")),
      thePt_Min(par.getParameter<double>("Pt_Min")),
      dzOptionString(par.getParameter<std::string>("dzOption")) {
  if (!dzOptionString.compare("dz"))
    dzOption = EgammaTrackSelector::dz;
  else if (!dzOptionString.compare("vz"))
    dzOption = EgammaTrackSelector::vz;
  else if (!dzOptionString.compare("bs"))
    dzOption = EgammaTrackSelector::bs;
  else if (!dzOptionString.compare("vtx"))
    dzOption = EgammaTrackSelector::vtx;
  else
    dzOption = EgammaTrackSelector::dz;
}

reco::IsoDeposit::Vetos EgammaTrackExtractor::vetos(const edm::Event& ev,
                                                    const edm::EventSetup& evSetup,
                                                    const reco::Track& track) const {
  reco::isodeposit::Direction dir(track.eta(), track.phi());
  return reco::IsoDeposit::Vetos(1, veto(dir));
}

reco::IsoDeposit::Veto EgammaTrackExtractor::veto(const reco::IsoDeposit::Direction& dir) const {
  reco::IsoDeposit::Veto result;
  result.vetoDir = dir;
  result.dR = theDR_Veto;
  return result;
}

IsoDeposit EgammaTrackExtractor::deposit(const Event& event,
                                         const EventSetup& eventSetup,
                                         const Candidate& candTk) const {
  static const std::string metname = "EgammaIsolationAlgos|EgammaTrackExtractor";

  reco::isodeposit::Direction candDir;
  double dzCut = 0;

  reco::TrackBase::Point beamPoint(0, 0, 0);
  if (theBeamlineOption == "BeamSpotFromEvent") {
    //pick beamSpot
    reco::BeamSpot beamSpot;
    edm::Handle<reco::BeamSpot> beamSpotH;

    event.getByToken(theBeamSpotToken, beamSpotH);

    if (beamSpotH.isValid()) {
      beamPoint = beamSpotH->position();
    }
  }

  Handle<View<Track> > tracksH;
  event.getByToken(theTrackCollectionToken, tracksH);

  if (candTk.isElectron()) {
    const reco::GsfElectron* elec = dynamic_cast<const reco::GsfElectron*>(&candTk);
    candDir = reco::isodeposit::Direction(elec->gsfTrack()->eta(), elec->gsfTrack()->phi());
  } else {
    candDir = reco::isodeposit::Direction(candTk.eta(), candTk.phi());
  }

  IsoDeposit deposit(candDir);
  deposit.setVeto(veto(candDir));
  deposit.addCandEnergy(candTk.et());

  View<Track>::const_iterator itrTr = tracksH->begin();
  View<Track>::const_iterator trEnd = tracksH->end();
  for (itrTr = tracksH->begin(); itrTr != trEnd; ++itrTr) {
    if (candDir.deltaR(reco::isodeposit::Direction(itrTr->eta(), itrTr->phi())) > theDR_Max)
      continue;

    if (itrTr->normalizedChi2() > theChi2Ndof_Max)
      continue;

    if (itrTr->pt() < thePt_Min)
      continue;

    if (theChi2Prob_Min > 0 && ChiSquaredProbability(itrTr->chi2(), itrTr->ndof()) < theChi2Prob_Min)
      continue;

    if (theNHits_Min > 0 && itrTr->numberOfValidHits() < theNHits_Min)
      continue;

    if (candTk.isElectron()) {
      const reco::GsfElectron* elec = dynamic_cast<const reco::GsfElectron*>(&candTk);
      switch (dzOption) {
        case EgammaTrackSelector::dz:
          dzCut = elec->gsfTrack()->dz() - itrTr->dz();
          break;
        case EgammaTrackSelector::vz:
          dzCut = elec->gsfTrack()->vz() - itrTr->vz();
          break;
        case EgammaTrackSelector::bs:
          dzCut = elec->gsfTrack()->dz(beamPoint) - itrTr->dz(beamPoint);
          break;
        case EgammaTrackSelector::vtx:
          dzCut = itrTr->dz(elec->gsfTrack()->vertex());
          break;
        default:
          dzCut = elec->gsfTrack()->vz() - itrTr->vz();
          break;
      }
    } else {
      switch (dzOption) {
        case EgammaTrackSelector::dz:
          dzCut = (*itrTr).dz() - candTk.vertex().z();
          break;
        case EgammaTrackSelector::vz:
          dzCut = (*itrTr).vz() - candTk.vertex().z();
          break;
        case EgammaTrackSelector::bs:
          dzCut = (*itrTr).dz(beamPoint) - candTk.vertex().z();
          break;
        case EgammaTrackSelector::vtx:
          dzCut = (*itrTr).dz(candTk.vertex());
          break;
        default:
          dzCut = (*itrTr).vz() - candTk.vertex().z();
          break;
      }
    }

    if (fabs(dzCut) > theDiff_z)
      continue;

    if (fabs(itrTr->dxy(beamPoint)) > theDiff_r)
      continue;

    deposit.addDeposit(reco::isodeposit::Direction(itrTr->eta(), itrTr->phi()), itrTr->pt());
  }

  return deposit;
}