CutBasedElectronID.cc

Go to the documentation of this file.

#include "RecoEgamma/ElectronIdentification/interface/CutBasedElectronID.h"
#include "DataFormats/EgammaReco/interface/BasicCluster.h"
#include "DataFormats/GsfTrackReco/interface/GsfTrack.h"
#include "DataFormats/BeamSpot/interface/BeamSpot.h"

#include <algorithm>

CutBasedElectronID::CutBasedElectronID(const edm::ParameterSet& conf, edm::ConsumesCollector& iC) {
  verticesCollection_ = iC.consumes<std::vector<reco::Vertex> >(conf.getParameter<edm::InputTag>("verticesCollection"));
}

void CutBasedElectronID::setup(const edm::ParameterSet& conf) {
  // Get all the parameters
  //baseSetup(conf);

  type_ = conf.getParameter<std::string>("electronIDType");
  quality_ = conf.getParameter<std::string>("electronQuality");
  version_ = conf.getParameter<std::string>("electronVersion");
  //verticesCollection_ = conf.getParameter<edm::InputTag>("verticesCollection");

  if (type_ == "classbased" and (version_ == "V06")) {
    newCategories_ = conf.getParameter<bool>("additionalCategories");
  }

  if (type_ == "classbased" and (version_ == "V03" or version_ == "V04" or version_ == "V05" or version_.empty())) {
    wantBinning_ = conf.getParameter<bool>("etBinning");
    newCategories_ = conf.getParameter<bool>("additionalCategories");
  }

  if (type_ == "robust" || type_ == "classbased") {
    std::string stringCut = type_ + quality_ + "EleIDCuts" + version_;
    cuts_ = conf.getParameter<edm::ParameterSet>(stringCut);
  } else {
    throw cms::Exception("Configuration")
        << "Invalid electronType parameter in CutBasedElectronID: must be robust or classbased\n";
  }
}

double CutBasedElectronID::result(const reco::GsfElectron* electron, const edm::Event& e, const edm::EventSetup& es) {
  if (type_ == "classbased")
    return cicSelection(electron, e, es);
  else if (type_ == "robust")
    return robustSelection(electron, e, es);

  return 0;
}

int CutBasedElectronID::classify(const reco::GsfElectron* electron) {
  double eta = fabs(electron->superCluster()->eta());
  double eOverP = electron->eSuperClusterOverP();
  double fBrem = electron->fbrem();

  int cat = -1;
  if (version_ == "V00" || version_ == "V01") {
    if ((electron->isEB() && fBrem < 0.06) || (electron->isEE() && fBrem < 0.1))
      cat = 1;
    else if (eOverP < 1.2 && eOverP > 0.8)
      cat = 0;
    else
      cat = 2;

    return cat;

  } else if (version_ == "V02") {
    if (electron->isEB()) {  // BARREL
      if (fBrem < 0.12)
        cat = 1;
      else if (eOverP < 1.2 && eOverP > 0.9)
        cat = 0;
      else
        cat = 2;
    } else {  // ENDCAP
      if (fBrem < 0.2)
        cat = 1;
      else if (eOverP < 1.22 && eOverP > 0.82)
        cat = 0;
      else
        cat = 2;
    }

    return cat;

  } else {
    if (electron->isEB()) {
      if ((fBrem >= 0.12) and (eOverP > 0.9) and (eOverP < 1.2))
        cat = 0;
      else if (((eta > .445 and eta < .45) or (eta > .79 and eta < .81) or (eta > 1.137 and eta < 1.157) or
                (eta > 1.47285 and eta < 1.4744)) and
               newCategories_)
        cat = 6;
      else if (electron->trackerDrivenSeed() and !electron->ecalDrivenSeed() and newCategories_)
        cat = 8;
      else if (fBrem < 0.12)
        cat = 1;
      else
        cat = 2;
    } else {
      if ((fBrem >= 0.2) and (eOverP > 0.82) and (eOverP < 1.22))
        cat = 3;
      else if (eta > 1.5 and eta < 1.58 and newCategories_)
        cat = 7;
      else if (electron->trackerDrivenSeed() and !electron->ecalDrivenSeed() and newCategories_)
        cat = 8;
      else if (fBrem < 0.2)
        cat = 4;
      else
        cat = 5;
    }

    return cat;
  }

  return -1;
}

double CutBasedElectronID::cicSelection(const reco::GsfElectron* electron,
                                        const edm::Event& e,
                                        const edm::EventSetup& es) {
  double scTheta = (2 * atan(exp(-electron->superCluster()->eta())));
  double scEt = electron->superCluster()->energy() * sin(scTheta);

  double eta = fabs(electron->superCluster()->eta());

  double eOverP = electron->eSuperClusterOverP();
  double eSeedOverPin = electron->eSeedClusterOverP();
  double fBrem = electron->fbrem();
  double hOverE = electron->hadronicOverEm();
  double sigmaee = electron->sigmaIetaIeta();  //sqrt(vLocCov[0]);
  double deltaPhiIn = electron->deltaPhiSuperClusterTrackAtVtx();
  double deltaEtaIn = electron->deltaEtaSuperClusterTrackAtVtx();

  double ip = 0;
  int mishits = electron->gsfTrack()->hitPattern().numberOfLostHits(reco::HitPattern::MISSING_INNER_HITS);
  double tkIso = electron->dr03TkSumPt();
  double ecalIso = electron->dr04EcalRecHitSumEt();
  double hcalIso = electron->dr04HcalTowerSumEt();

  if (version_ == "V00") {
    sigmaee = electron->sigmaEtaEta();  //sqrt(vCov[0]);
    if (electron->isEE())
      sigmaee = sigmaee - 0.02 * (fabs(eta) - 2.3);  //correct sigmaetaeta dependence on eta in endcap
  }

  if (version_ != "V01" or version_ != "V00") {
    edm::Handle<reco::VertexCollection> vtxH;
    e.getByToken(verticesCollection_, vtxH);
    if (!vtxH->empty()) {
      reco::VertexRef vtx(vtxH, 0);
      ip = fabs(electron->gsfTrack()->dxy(math::XYZPoint(vtx->x(), vtx->y(), vtx->z())));
    } else
      ip = fabs(electron->gsfTrack()->dxy());

    if (electron->isEB()) {
      sigmaee = electron->sigmaIetaIeta();  //sqrt(vCov[0]);
    }
  }

  std::vector<double> cut;

  int cat = classify(electron);
  int eb;

  if (electron->isEB())
    eb = 0;
  else
    eb = 1;

  // LOOSE and TIGHT Selections
  if (type_ == "classbased" && (version_ == "V01" || version_ == "V00")) {
    if ((eOverP < 0.8) && (fBrem < 0.2))
      return 0.;

    cut = cuts_.getParameter<std::vector<double> >("hOverE");
    if (hOverE > cut[cat + 4 * eb])
      return 0.;

    cut = cuts_.getParameter<std::vector<double> >("sigmaEtaEta");
    if (sigmaee > cut[cat + 4 * eb])
      return 0.;

    cut = cuts_.getParameter<std::vector<double> >("deltaPhiIn");
    if (eOverP < 1.5) {
      if (fabs(deltaPhiIn) > cut[cat + 4 * eb])
        return 0.;
    } else {
      if (fabs(deltaPhiIn) > cut[3 + 4 * eb])
        return 0.;
    }

    cut = cuts_.getParameter<std::vector<double> >("deltaEtaIn");
    if (fabs(deltaEtaIn) > cut[cat + 4 * eb])
      return 0.;

    cut = cuts_.getParameter<std::vector<double> >("eSeedOverPin");
    if (eSeedOverPin < cut[cat + 4 * eb])
      return 0.;

    if (quality_ == "tight")
      if (eOverP < 0.9 * (1 - fBrem))
        return 0.;

    return 1.;
  }

  if (type_ == "classbased" and version_ == "V02") {
    double result = 0.;

    int bin = 0;

    if (scEt < 20.)
      bin = 2;
    else if (scEt > 30.)
      bin = 0;
    else
      bin = 1;

    if (fBrem > 0)
      eSeedOverPin = eSeedOverPin + fBrem;

    if (bin != 2) {
      tkIso = tkIso * pow(40. / scEt, 2);
      ecalIso = ecalIso * pow(40. / scEt, 2);
      hcalIso = hcalIso * pow(40. / scEt, 2);
    }

    std::vector<double> cutTk = cuts_.getParameter<std::vector<double> >("cutisotk");
    std::vector<double> cutEcal = cuts_.getParameter<std::vector<double> >("cutisoecal");
    std::vector<double> cutHcal = cuts_.getParameter<std::vector<double> >("cutisohcal");
    if ((tkIso > cutTk[cat + 3 * eb + bin * 6]) || (ecalIso > cutEcal[cat + 3 * eb + bin * 6]) ||
        (hcalIso > cutHcal[cat + 3 * eb + bin * 6]))
      result = 0.;
    else
      result = 2.;

    if (fBrem > -2) {
      std::vector<double> cuthoe = cuts_.getParameter<std::vector<double> >("cuthoe");
      std::vector<double> cutsee = cuts_.getParameter<std::vector<double> >("cutsee");
      std::vector<double> cutdphi = cuts_.getParameter<std::vector<double> >("cutdphiin");
      std::vector<double> cutdeta = cuts_.getParameter<std::vector<double> >("cutdetain");
      std::vector<double> cuteopin = cuts_.getParameter<std::vector<double> >("cuteseedopcor");
      std::vector<double> cutet = cuts_.getParameter<std::vector<double> >("cutet");
      std::vector<double> cutip = cuts_.getParameter<std::vector<double> >("cutip");
      std::vector<double> cutmishits = cuts_.getParameter<std::vector<double> >("cutmishits");
      if ((hOverE < cuthoe[cat + 3 * eb + bin * 6]) and (sigmaee < cutsee[cat + 3 * eb + bin * 6]) and
          (fabs(deltaPhiIn) < cutdphi[cat + 3 * eb + bin * 6]) and
          (fabs(deltaEtaIn) < cutdeta[cat + 3 * eb + bin * 6]) and (eSeedOverPin > cuteopin[cat + 3 * eb + bin * 6]) and
          (ip < cutip[cat + 3 * eb + bin * 6]) and (mishits < cutmishits[cat + 3 * eb + bin * 6]))
        result = result + 1.;
    }
    return result;
  }

  if (version_ == "V03" or version_ == "V04" or version_ == "V05") {
    double result = 0.;

    int bin = 0;

    if (wantBinning_) {
      if (scEt < 20.)
        bin = 2;
      else if (scEt > 30.)
        bin = 0;
      else
        bin = 1;
    }

    if (fBrem > 0)
      eSeedOverPin = eSeedOverPin + fBrem;

    float iso_sum = tkIso + ecalIso + hcalIso;
    float iso_sum_corrected = iso_sum * pow(40. / scEt, 2);

    std::vector<double> cutIsoSum = cuts_.getParameter<std::vector<double> >("cutiso_sum");
    std::vector<double> cutIsoSumCorr = cuts_.getParameter<std::vector<double> >("cutiso_sumoet");
    if ((iso_sum < cutIsoSum[cat + bin * 9]) and (iso_sum_corrected < cutIsoSumCorr[cat + bin * 9]))
      result += 2.;

    if (fBrem > -2) {
      std::vector<double> cuthoe = cuts_.getParameter<std::vector<double> >("cuthoe");
      std::vector<double> cutsee = cuts_.getParameter<std::vector<double> >("cutsee");
      std::vector<double> cutdphi = cuts_.getParameter<std::vector<double> >("cutdphiin");
      std::vector<double> cutdeta = cuts_.getParameter<std::vector<double> >("cutdetain");
      std::vector<double> cuteopin = cuts_.getParameter<std::vector<double> >("cuteseedopcor");
      std::vector<double> cutet = cuts_.getParameter<std::vector<double> >("cutet");

      if ((hOverE < cuthoe[cat + bin * 9]) and (sigmaee < cutsee[cat + bin * 9]) and
          (fabs(deltaPhiIn) < cutdphi[cat + bin * 9]) and (fabs(deltaEtaIn) < cutdeta[cat + bin * 9]) and
          (eSeedOverPin > cuteopin[cat + bin * 9]) and (scEt > cutet[cat + bin * 9]))
        result += 1.;
    }

    std::vector<double> cutip = cuts_.getParameter<std::vector<double> >("cutip_gsf");
    if (ip < cutip[cat + bin * 9])
      result += 8;

    std::vector<double> cutmishits = cuts_.getParameter<std::vector<double> >("cutfmishits");
    std::vector<double> cutdcotdist = cuts_.getParameter<std::vector<double> >("cutdcotdist");

    float dist = (electron->convDist() == -9999. ? 9999 : electron->convDist());
    float dcot = (electron->convDcot() == -9999. ? 9999 : electron->convDcot());

    float dcotdistcomb =
        ((0.04 - std::max(fabs(dist), fabs(dcot))) > 0 ? (0.04 - std::max(fabs(dist), fabs(dcot))) : 0);

    if ((mishits < cutmishits[cat + bin * 9]) and (dcotdistcomb < cutdcotdist[cat + bin * 9]))
      result += 4;

    return result;
  }

  if (type_ == "classbased" && (version_ == "V06" || version_.empty())) {
    std::vector<double> cutIsoSum = cuts_.getParameter<std::vector<double> >("cutiso_sum");
    std::vector<double> cutIsoSumCorr = cuts_.getParameter<std::vector<double> >("cutiso_sumoet");
    std::vector<double> cuthoe = cuts_.getParameter<std::vector<double> >("cuthoe");
    std::vector<double> cutsee = cuts_.getParameter<std::vector<double> >("cutsee");
    std::vector<double> cutdphi = cuts_.getParameter<std::vector<double> >("cutdphiin");
    std::vector<double> cutdeta = cuts_.getParameter<std::vector<double> >("cutdetain");
    std::vector<double> cuteopin = cuts_.getParameter<std::vector<double> >("cuteseedopcor");
    std::vector<double> cutmishits = cuts_.getParameter<std::vector<double> >("cutfmishits");
    std::vector<double> cutdcotdist = cuts_.getParameter<std::vector<double> >("cutdcotdist");
    std::vector<double> cutip = cuts_.getParameter<std::vector<double> >("cutip_gsf");
    std::vector<double> cutIsoSumCorrl = cuts_.getParameter<std::vector<double> >("cutiso_sumoetl");
    std::vector<double> cuthoel = cuts_.getParameter<std::vector<double> >("cuthoel");
    std::vector<double> cutseel = cuts_.getParameter<std::vector<double> >("cutseel");
    std::vector<double> cutdphil = cuts_.getParameter<std::vector<double> >("cutdphiinl");
    std::vector<double> cutdetal = cuts_.getParameter<std::vector<double> >("cutdetainl");
    std::vector<double> cutipl = cuts_.getParameter<std::vector<double> >("cutip_gsfl");

    int result = 0;

    const int ncuts = 10;
    std::vector<bool> cut_results(ncuts, false);

    float iso_sum = tkIso + ecalIso + hcalIso;
    float scEta = electron->superCluster()->eta();
    if (fabs(scEta) > 1.5)
      iso_sum += (fabs(scEta) - 1.5) * 1.09;

    float iso_sumoet = iso_sum * (40. / scEt);

    float eseedopincor = eSeedOverPin + fBrem;
    if (fBrem < 0)
      eseedopincor = eSeedOverPin;

    float dist = (electron->convDist() == -9999. ? 9999 : electron->convDist());
    float dcot = (electron->convDcot() == -9999. ? 9999 : electron->convDcot());

    float dcotdistcomb =
        ((0.04 - std::max(fabs(dist), fabs(dcot))) > 0 ? (0.04 - std::max(fabs(dist), fabs(dcot))) : 0);

    for (int cut = 0; cut < ncuts; cut++) {
      switch (cut) {
        case 0:
          cut_results[cut] = compute_cut(fabs(deltaEtaIn), scEt, cutdetal[cat], cutdeta[cat]);
          break;
        case 1:
          cut_results[cut] = compute_cut(fabs(deltaPhiIn), scEt, cutdphil[cat], cutdphi[cat]);
          break;
        case 2:
          cut_results[cut] = (eseedopincor > cuteopin[cat]);
          break;
        case 3:
          cut_results[cut] = compute_cut(hOverE, scEt, cuthoel[cat], cuthoe[cat]);
          break;
        case 4:
          cut_results[cut] = compute_cut(sigmaee, scEt, cutseel[cat], cutsee[cat]);
          break;
        case 5:
          cut_results[cut] = compute_cut(iso_sumoet, scEt, cutIsoSumCorrl[cat], cutIsoSumCorr[cat]);
          break;
        case 6:
          cut_results[cut] = (iso_sum < cutIsoSum[cat]);
          break;
        case 7:
          cut_results[cut] = compute_cut(fabs(ip), scEt, cutipl[cat], cutip[cat]);
          break;
        case 8:
          cut_results[cut] = (mishits < cutmishits[cat]);
          break;
        case 9:
          cut_results[cut] = (dcotdistcomb < cutdcotdist[cat]);
          break;
      }
    }

    // ID part
    if (cut_results[0] & cut_results[1] & cut_results[2] & cut_results[3] & cut_results[4])
      result = result + 1;

    // ISO part
    if (cut_results[5] & cut_results[6])
      result = result + 2;

    // IP part
    if (cut_results[7])
      result = result + 8;

    // Conversion part
    if (cut_results[8] & cut_results[9])
      result = result + 4;

    return result;
  }

  return -1.;
}

bool CutBasedElectronID::compute_cut(double x, double et, double cut_min, double cut_max, bool gtn) {
  float et_min = 10;
  float et_max = 40;

  bool accept = false;
  float cut = cut_max;  //  the cut at et=40 GeV

  if (et < et_max) {
    cut = cut_min + (1 / et_min - 1 / et) * (cut_max - cut_min) / (1 / et_min - 1 / et_max);
  }

  if (et < et_min) {
    cut = cut_min;
  }

  if (gtn) {  // useful for e/p cut which is gt
    accept = (x >= cut);
  } else {
    accept = (x <= cut);
  }

  //std::cout << x << " " << cut_min << " " << cut << " " << cut_max << " " << et << " " << accept << std::endl;
  return accept;
}

double CutBasedElectronID::robustSelection(const reco::GsfElectron* electron,
                                           const edm::Event& e,
                                           const edm::EventSetup& es) {
  double scTheta = (2 * atan(exp(-electron->superCluster()->eta())));
  double scEt = electron->superCluster()->energy() * sin(scTheta);
  double eta = electron->p4().Eta();
  double eOverP = electron->eSuperClusterOverP();
  double hOverE = electron->hadronicOverEm();
  double sigmaee = electron->sigmaIetaIeta();
  double e25Max = electron->e2x5Max();
  double e15 = electron->e1x5();
  double e55 = electron->e5x5();
  double e25Maxoe55 = e25Max / e55;
  double e15oe55 = e15 / e55;
  double deltaPhiIn = electron->deltaPhiSuperClusterTrackAtVtx();
  double deltaEtaIn = electron->deltaEtaSuperClusterTrackAtVtx();

  double ip = 0;
  int mishits = electron->gsfTrack()->hitPattern().numberOfLostHits(reco::HitPattern::MISSING_INNER_HITS);
  double tkIso = electron->dr03TkSumPt();
  double ecalIso = electron->dr04EcalRecHitSumEt();
  double ecalIsoPed = (electron->isEB()) ? std::max(0., ecalIso - 1.) : ecalIso;
  double hcalIso = electron->dr04HcalTowerSumEt();
  double hcalIso1 = electron->dr04HcalDepth1TowerSumEt();
  double hcalIso2 = electron->dr04HcalDepth2TowerSumEt();

  if (version_ == "V00") {
    sigmaee = electron->sigmaEtaEta();
    if (electron->isEE())
      sigmaee = sigmaee - 0.02 * (fabs(eta) - 2.3);  //correct sigmaetaeta dependence on eta in endcap
  }

  if (version_ == "V03" or version_ == "V04") {
    edm::Handle<reco::BeamSpot> pBeamSpot;
    // uses the same name for the vertex collection to avoid adding more new names
    e.getByToken(verticesCollection_, pBeamSpot);
    if (pBeamSpot.isValid()) {
      const reco::BeamSpot* bspot = pBeamSpot.product();
      const math::XYZPoint& bspotPosition = bspot->position();
      ip = fabs(electron->gsfTrack()->dxy(bspotPosition));
    } else
      ip = fabs(electron->gsfTrack()->dxy());
  }

  if (version_ == "V04" or version_ == "V05") {
    ecalIso = electron->dr03EcalRecHitSumEt();
    ecalIsoPed = (electron->isEB()) ? std::max(0., ecalIso - 1.) : ecalIso;
    hcalIso = electron->dr03HcalTowerSumEt();
    hcalIso1 = electron->dr03HcalDepth1TowerSumEt();
    hcalIso2 = electron->dr03HcalDepth2TowerSumEt();
  }

  if (version_ == "V05") {
    edm::Handle<reco::VertexCollection> vtxH;
    e.getByToken(verticesCollection_, vtxH);
    if (!vtxH->empty()) {
      reco::VertexRef vtx(vtxH, 0);
      ip = fabs(electron->gsfTrack()->dxy(math::XYZPoint(vtx->x(), vtx->y(), vtx->z())));
    } else
      ip = fabs(electron->gsfTrack()->dxy());
  }

  // .....................................................................................
  std::vector<double> cut;
  // ROBUST Selection
  if (type_ == "robust") {
    double result = 0;

    // hoe, sigmaEtaEta, dPhiIn, dEtaIn
    if (electron->isEB())
      cut = cuts_.getParameter<std::vector<double> >("barrel");
    else
      cut = cuts_.getParameter<std::vector<double> >("endcap");
    // check isolations: if only isolation passes result = 2
    if (quality_ == "highenergy") {
      if ((tkIso > cut[6] || hcalIso2 > cut[11]) ||
          (electron->isEB() && ((ecalIso + hcalIso1) > cut[7] + cut[8] * scEt)) ||
          (electron->isEE() && (scEt >= 50.) && ((ecalIso + hcalIso1) > cut[7] + cut[8] * (scEt - 50))) ||
          (electron->isEE() && (scEt < 50.) && ((ecalIso + hcalIso1) > cut[9] + cut[10] * (scEt - 50))))
        result = 0;
      else
        result = 2;
    } else {
      if ((tkIso > cut[6]) || (ecalIso > cut[7]) || (hcalIso > cut[8]) || (hcalIso1 > cut[9]) || (hcalIso2 > cut[10]) ||
          (tkIso / electron->p4().Pt() > cut[11]) || (ecalIso / electron->p4().Pt() > cut[12]) ||
          (hcalIso / electron->p4().Pt() > cut[13]) || ((tkIso + ecalIso + hcalIso) > cut[14]) ||
          (((tkIso + ecalIso + hcalIso) / electron->p4().Pt()) > cut[15]) ||
          ((tkIso + ecalIsoPed + hcalIso) > cut[16]) ||
          (((tkIso + ecalIsoPed + hcalIso) / electron->p4().Pt()) > cut[17]))
        result = 0.;
      else
        result = 2.;
    }

    if ((hOverE < cut[0]) && (sigmaee < cut[1]) && (fabs(deltaPhiIn) < cut[2]) && (fabs(deltaEtaIn) < cut[3]) &&
        (e25Maxoe55 > cut[4] && e15oe55 > cut[5]) && (sigmaee >= cut[18]) && (eOverP > cut[19] && eOverP < cut[20])) {
      result = result + 1;
    }

    if (ip > cut[21])
      return result;
    if (mishits > cut[22])  // expected missing hits
      return result;
    // positive cut[23] means to demand a valid hit in 1st layer PXB
    if (cut[23] > 0 &&
        !electron->gsfTrack()->hitPattern().hasValidHitInPixelLayer(PixelSubdetector::SubDetector::PixelBarrel, 1))
      return result;

    // cut[24]: Dist cut[25]: dcot
    float dist = fabs(electron->convDist());
    float dcot = fabs(electron->convDcot());
    bool isConversion = (cut[24] > 99. || cut[25] > 99.) ? false : (dist < cut[24] && dcot < cut[25]);
    if (isConversion)
      return result;

    result += 4;

    return result;
  }

  return -1.;
}