EleTkIsolFromCands.cc

Go to the documentation of this file.

#include "RecoEgamma/EgammaIsolationAlgos/interface/EleTkIsolFromCands.h"
#include "DataFormats/TrackReco/interface/Track.h"
#include "DataFormats/Math/interface/deltaR.h"

EleTkIsolFromCands::TrkCuts::TrkCuts(const edm::ParameterSet& para) {
  minPt = para.getParameter<double>("minPt");
  auto sq = [](double val) { return val * val; };
  minDR2 = sq(para.getParameter<double>("minDR"));
  maxDR2 = sq(para.getParameter<double>("maxDR"));
  minDEta = para.getParameter<double>("minDEta");
  maxDZ = para.getParameter<double>("maxDZ");
  minHits = para.getParameter<int>("minHits");
  minPixelHits = para.getParameter<int>("minPixelHits");
  maxDPtPt = para.getParameter<double>("maxDPtPt");

  auto qualNames = para.getParameter<std::vector<std::string> >("allowedQualities");
  auto algoNames = para.getParameter<std::vector<std::string> >("algosToReject");

  for (auto& qualName : qualNames) {
    allowedQualities.push_back(reco::TrackBase::qualityByName(qualName));
  }
  for (auto& algoName : algoNames) {
    algosToReject.push_back(reco::TrackBase::algoByName(algoName));
  }
  std::sort(algosToReject.begin(), algosToReject.end());
}

edm::ParameterSetDescription EleTkIsolFromCands::TrkCuts::pSetDescript() {
  edm::ParameterSetDescription desc;
  desc.add<double>("minPt", 1.0);
  desc.add<double>("maxDR", 0.3);
  desc.add<double>("minDR", 0.000);
  desc.add<double>("minDEta", 0.005);
  desc.add<double>("maxDZ", 0.1);
  desc.add<double>("maxDPtPt", -1);
  desc.add<int>("minHits", 8);
  desc.add<int>("minPixelHits", 1);
  desc.add<std::vector<std::string> >("allowedQualities");
  desc.add<std::vector<std::string> >("algosToReject");
  return desc;
}

EleTkIsolFromCands::EleTkIsolFromCands(const edm::ParameterSet& para)
    : barrelCuts_(para.getParameter<edm::ParameterSet>("barrelCuts")),
      endcapCuts_(para.getParameter<edm::ParameterSet>("endcapCuts")) {}

edm::ParameterSetDescription EleTkIsolFromCands::pSetDescript() {
  edm::ParameterSetDescription desc;
  desc.add("barrelCuts", TrkCuts::pSetDescript());
  desc.add("endcapCuts", TrkCuts::pSetDescript());
  return desc;
}

std::pair<int, double> EleTkIsolFromCands::calIsol(const reco::TrackBase& eleTrk,
                                                   const pat::PackedCandidateCollection& cands,
                                                   const PIDVeto pidVeto) const {
  return calIsol(eleTrk.eta(), eleTrk.phi(), eleTrk.vz(), cands, pidVeto);
}

std::pair<int, double> EleTkIsolFromCands::calIsol(const double eleEta,
                                                   const double elePhi,
                                                   const double eleVZ,
                                                   const pat::PackedCandidateCollection& cands,
                                                   const PIDVeto pidVeto) const {
  double ptSum = 0.;
  int nrTrks = 0;

  const TrkCuts& cuts = std::abs(eleEta) < 1.5 ? barrelCuts_ : endcapCuts_;

  for (auto& cand : cands) {
    if (cand.hasTrackDetails() && cand.charge() != 0 && passPIDVeto(cand.pdgId(), pidVeto)) {
      const reco::Track& trk = cand.pseudoTrack();
      if (passTrkSel(trk, trk.pt(), cuts, eleEta, elePhi, eleVZ)) {
        ptSum += trk.pt();
        nrTrks++;
      }
    }
  }
  return {nrTrks, ptSum};
}

std::pair<int, double> EleTkIsolFromCands::calIsol(const reco::TrackBase& eleTrk,
                                                   const reco::TrackCollection& tracks) const {
  return calIsol(eleTrk.eta(), eleTrk.phi(), eleTrk.vz(), tracks);
}

std::pair<int, double> EleTkIsolFromCands::calIsol(const double eleEta,
                                                   const double elePhi,
                                                   const double eleVZ,
                                                   const reco::TrackCollection& tracks) const {
  double ptSum = 0.;
  int nrTrks = 0;

  const TrkCuts& cuts = std::abs(eleEta) < 1.5 ? barrelCuts_ : endcapCuts_;

  for (auto& trk : tracks) {
    if (passTrkSel(trk, trk.pt(), cuts, eleEta, elePhi, eleVZ)) {
      ptSum += trk.pt();
      nrTrks++;
    }
  }
  return {nrTrks, ptSum};
}

bool EleTkIsolFromCands::passPIDVeto(const int pdgId, const EleTkIsolFromCands::PIDVeto veto) {
  int pidAbs = std::abs(pdgId);
  switch (veto) {
    case PIDVeto::NONE:
      return true;
    case PIDVeto::ELES:
      if (pidAbs == 11)
        return false;
      else
        return true;
    case PIDVeto::NONELES:
      if (pidAbs == 11)
        return true;
      else
        return false;
  }
  throw cms::Exception("CodeError") << "invalid PIDVeto " << static_cast<int>(veto) << ", "
                                    << "this is likely due to some static casting of invalid ints somewhere";
}

EleTkIsolFromCands::PIDVeto EleTkIsolFromCands::pidVetoFromStr(const std::string& vetoStr) {
  if (vetoStr == "NONE")
    return PIDVeto::NONE;
  else if (vetoStr == "ELES")
    return PIDVeto::ELES;
  else if (vetoStr == "NONELES")
    return PIDVeto::NONELES;
  else {
    throw cms::Exception("CodeError") << "unrecognised string " << vetoStr
                                      << ", either a typo or this function needs to be updated";
  }
}

bool EleTkIsolFromCands::passTrkSel(const reco::TrackBase& trk,
                                    const double trkPt,
                                    const TrkCuts& cuts,
                                    const double eleEta,
                                    const double elePhi,
                                    const double eleVZ) {
  const float dR2 = reco::deltaR2(eleEta, elePhi, trk.eta(), trk.phi());
  const float dEta = trk.eta() - eleEta;
  const float dZ = eleVZ - trk.vz();

  return dR2 >= cuts.minDR2 && dR2 <= cuts.maxDR2 && std::abs(dEta) >= cuts.minDEta && std::abs(dZ) < cuts.maxDZ &&
         trk.hitPattern().numberOfValidHits() >= cuts.minHits &&
         trk.hitPattern().numberOfValidPixelHits() >= cuts.minPixelHits &&
         (trk.ptError() / trkPt < cuts.maxDPtPt || cuts.maxDPtPt < 0) && passQual(trk, cuts.allowedQualities) &&
         passAlgo(trk, cuts.algosToReject) && trkPt > cuts.minPt;
}

bool EleTkIsolFromCands::passQual(const reco::TrackBase& trk, const std::vector<reco::TrackBase::TrackQuality>& quals) {
  if (quals.empty())
    return true;

  for (auto qual : quals) {
    if (trk.quality(qual))
      return true;
  }

  return false;
}

bool EleTkIsolFromCands::passAlgo(const reco::TrackBase& trk,
                                  const std::vector<reco::TrackBase::TrackAlgorithm>& algosToRej) {
  return algosToRej.empty() || !std::binary_search(algosToRej.begin(), algosToRej.end(), trk.algo());
}