MvaMEtUtilities.cc

Go to the documentation of this file.

#include "RecoMET/METPUSubtraction/interface/MvaMEtUtilities.h"

#include "DataFormats/Math/interface/deltaR.h"

#include <algorithm>
#include <cmath>

MvaMEtUtilities::MvaMEtUtilities(const edm::ParameterSet& cfg) {
  // jet id
  /* ===> this code parses an xml it uses parameter set 
  edm::ParameterSet jetConfig = iConfig.getParameter<edm::ParameterSet>("JetIdParams");
  for(int i0 = 0; i0 < 3; i0++) { 
    std::string lCutType                            = "Tight";
    if(i0 == PileupJetIdentifier::kMedium) lCutType = "Medium";
    if(i0 == PileupJetIdentifier::kLoose)  lCutType = "Loose";
    std::vector<double> pt010  = jetConfig.getParameter<std::vector<double> >(("Pt010_" +lCutType).c_str());
    std::vector<double> pt1020 = jetConfig.getParameter<std::vector<double> >(("Pt1020_"+lCutType).c_str());
    std::vector<double> pt2030 = jetConfig.getParameter<std::vector<double> >(("Pt2030_"+lCutType).c_str());
    std::vector<double> pt3050 = jetConfig.getParameter<std::vector<double> >(("Pt3050_"+lCutType).c_str());
    for(int i1 = 0; i1 < 4; i1++) mvacut_[i0][0][i1] = pt010 [i1];
    for(int i1 = 0; i1 < 4; i1++) mvacut_[i0][1][i1] = pt1020[i1];
    for(int i1 = 0; i1 < 4; i1++) mvacut_[i0][2][i1] = pt2030[i1];
    for(int i1 = 0; i1 < 4; i1++) mvacut_[i0][3][i1] = pt3050[i1];
    }
  */
  //Tight Id => not used
  mvaCut_[0][0][0] = 0.5;
  mvaCut_[0][0][1] = 0.6;
  mvaCut_[0][0][2] = 0.6;
  mvaCut_[0][0][3] = 0.9;
  mvaCut_[0][1][0] = -0.2;
  mvaCut_[0][1][1] = 0.2;
  mvaCut_[0][1][2] = 0.2;
  mvaCut_[0][1][3] = 0.6;
  mvaCut_[0][2][0] = 0.3;
  mvaCut_[0][2][1] = 0.4;
  mvaCut_[0][2][2] = 0.7;
  mvaCut_[0][2][3] = 0.8;
  mvaCut_[0][3][0] = 0.5;
  mvaCut_[0][3][1] = 0.4;
  mvaCut_[0][3][2] = 0.8;
  mvaCut_[0][3][3] = 0.9;
  //Medium id => not used
  mvaCut_[1][0][0] = 0.2;
  mvaCut_[1][0][1] = 0.4;
  mvaCut_[1][0][2] = 0.2;
  mvaCut_[1][0][3] = 0.6;
  mvaCut_[1][1][0] = -0.3;
  mvaCut_[1][1][1] = 0.;
  mvaCut_[1][1][2] = 0.;
  mvaCut_[1][1][3] = 0.5;
  mvaCut_[1][2][0] = 0.2;
  mvaCut_[1][2][1] = 0.2;
  mvaCut_[1][2][2] = 0.5;
  mvaCut_[1][2][3] = 0.7;
  mvaCut_[1][3][0] = 0.3;
  mvaCut_[1][3][1] = 0.2;
  mvaCut_[1][3][2] = 0.7;
  mvaCut_[1][3][3] = 0.8;
  //Met Id => used
  mvaCut_[2][0][0] = -0.2;
  mvaCut_[2][0][1] = -0.3;
  mvaCut_[2][0][2] = -0.5;
  mvaCut_[2][0][3] = -0.5;
  mvaCut_[2][1][0] = -0.2;
  mvaCut_[2][1][1] = -0.2;
  mvaCut_[2][1][2] = -0.5;
  mvaCut_[2][1][3] = -0.3;
  mvaCut_[2][2][0] = -0.2;
  mvaCut_[2][2][1] = -0.2;
  mvaCut_[2][2][2] = -0.2;
  mvaCut_[2][2][3] = 0.1;
  mvaCut_[2][3][0] = -0.2;
  mvaCut_[2][3][1] = -0.2;
  mvaCut_[2][3][2] = 0.;
  mvaCut_[2][3][3] = 0.2;

  dzCut_ = cfg.getParameter<double>("dZcut");
  ptThreshold_ = (cfg.exists("ptThreshold")) ? cfg.getParameter<int>("ptThreshold") : -1000;
}

MvaMEtUtilities::~MvaMEtUtilities() {
  // nothing to be done yet...
}

bool MvaMEtUtilities::passesMVA(const reco::Candidate::LorentzVector& jetP4, double mvaJetId) {
  int ptBin = 0;
  if (jetP4.pt() >= 10. && jetP4.pt() < 20.)
    ptBin = 1;
  if (jetP4.pt() >= 20. && jetP4.pt() < 30.)
    ptBin = 2;
  if (jetP4.pt() >= 30.)
    ptBin = 3;

  int etaBin = 0;
  if (std::abs(jetP4.eta()) >= 2.5 && std::abs(jetP4.eta()) < 2.75)
    etaBin = 1;
  if (std::abs(jetP4.eta()) >= 2.75 && std::abs(jetP4.eta()) < 3.0)
    etaBin = 2;
  if (std::abs(jetP4.eta()) >= 3.0 && std::abs(jetP4.eta()) < 5.0)
    etaBin = 3;

  return (mvaJetId > mvaCut_[2][ptBin][etaBin]);
}

reco::Candidate::LorentzVector MvaMEtUtilities::leadJetP4(const std::vector<reco::PUSubMETCandInfo>& jets) {
  return jetP4(jets, 0);
}

reco::Candidate::LorentzVector MvaMEtUtilities::subleadJetP4(const std::vector<reco::PUSubMETCandInfo>& jets) {
  return jetP4(jets, 1);
}

reco::Candidate::LorentzVector MvaMEtUtilities::jetP4(const std::vector<reco::PUSubMETCandInfo>& jets, unsigned idx) {
  reco::Candidate::LorentzVector retVal(0., 0., 0., 0.);
  if (idx < jets.size()) {
    std::vector<reco::PUSubMETCandInfo> jets_sorted = jets;
    std::sort(jets_sorted.rbegin(), jets_sorted.rend());
    retVal = jets_sorted[idx].p4();
  }
  return retVal;
}
unsigned MvaMEtUtilities::numJetsAboveThreshold(const std::vector<reco::PUSubMETCandInfo>& jets, double ptThreshold) {
  unsigned retVal = 0;
  for (std::vector<reco::PUSubMETCandInfo>::const_iterator jet = jets.begin(); jet != jets.end(); ++jet) {
    if (jet->p4().pt() > ptThreshold)
      ++retVal;
  }
  return retVal;
}
std::vector<reco::PUSubMETCandInfo> MvaMEtUtilities::cleanJets(const std::vector<reco::PUSubMETCandInfo>& jets,
                                                               const std::vector<reco::PUSubMETCandInfo>& leptons,
                                                               double ptThreshold,
                                                               double dRmatch) {
  double dR2match = dRmatch * dRmatch;
  std::vector<reco::PUSubMETCandInfo> retVal;
  for (std::vector<reco::PUSubMETCandInfo>::const_iterator jet = jets.begin(); jet != jets.end(); ++jet) {
    bool isOverlap = false;
    for (std::vector<reco::PUSubMETCandInfo>::const_iterator lepton = leptons.begin(); lepton != leptons.end();
         ++lepton) {
      if (deltaR2(jet->p4(), lepton->p4()) < dR2match)
        isOverlap = true;
    }
    if (jet->p4().pt() > ptThreshold && !isOverlap)
      retVal.push_back(*jet);
  }
  return retVal;
}

std::vector<reco::PUSubMETCandInfo> MvaMEtUtilities::cleanPFCands(
    const std::vector<reco::PUSubMETCandInfo>& pfCandidates,
    const std::vector<reco::PUSubMETCandInfo>& leptons,
    double dRmatch,
    bool invert) {
  double dR2match = dRmatch * dRmatch;
  std::vector<reco::PUSubMETCandInfo> retVal;
  for (std::vector<reco::PUSubMETCandInfo>::const_iterator pfCandidate = pfCandidates.begin();
       pfCandidate != pfCandidates.end();
       ++pfCandidate) {
    bool isOverlap = false;
    for (std::vector<reco::PUSubMETCandInfo>::const_iterator lepton = leptons.begin(); lepton != leptons.end();
         ++lepton) {
      if (deltaR2(pfCandidate->p4(), lepton->p4()) < dR2match)
        isOverlap = true;
    }
    if ((!isOverlap && !invert) || (isOverlap && invert))
      retVal.push_back(*pfCandidate);
  }
  return retVal;
}
void MvaMEtUtilities::finalize(CommonMETData& metData) {
  metData.met = sqrt(metData.mex * metData.mex + metData.mey * metData.mey);
  metData.mez = 0.;
  metData.phi = atan2(metData.mey, metData.mex);
}

CommonMETData MvaMEtUtilities::computeCandSum(int compKey,
                                              double dZmax,
                                              int dZflag,
                                              bool iCharged,
                                              bool mvaPassFlag,
                                              const std::vector<reco::PUSubMETCandInfo>& objects) {
  CommonMETData retVal;
  retVal.mex = 0.;
  retVal.mey = 0.;
  retVal.sumet = 0.;

  for (std::vector<reco::PUSubMETCandInfo>::const_iterator object = objects.begin(); object != objects.end();
       ++object) {
    double pFrac = 1;

    //pf candidates
    // dZcut
    //   maximum distance within which tracks are
    //considered to be associated to hard scatter vertex
    // dZflag
    //   0 : select charged PFCandidates originating from hard scatter vertex
    //   1 : select charged PFCandidates originating from pile-up vertices
    //   2 : select all PFCandidates
    if (compKey == MvaMEtUtilities::kPFCands) {
      if (object->dZ() < 0. && dZflag != 2)
        continue;
      if (object->dZ() > dZmax && dZflag == 0)
        continue;
      if (object->dZ() < dZmax && dZflag == 1)
        continue;
    }

    //leptons
    if (compKey == MvaMEtUtilities::kLeptons) {
      if (iCharged)
        pFrac = object->chargedEnFrac();
    }

    //jets
    if (compKey == MvaMEtUtilities::kJets) {
      bool passesMVAjetId = passesMVA(object->p4(), object->mva());

      if (passesMVAjetId && !mvaPassFlag)
        continue;
      if (!passesMVAjetId && mvaPassFlag)
        continue;

      pFrac = 1 - object->chargedEnFrac();  //neutral energy fraction
    }

    retVal.mex += object->p4().px() * pFrac;
    retVal.mey += object->p4().py() * pFrac;
    retVal.sumet += object->p4().pt() * pFrac;
  }

  finalize(retVal);
  return retVal;
}

CommonMETData MvaMEtUtilities::computeRecoil(int metType) {
  CommonMETData retVal;

  if (metType == MvaMEtUtilities::kPF) {
    //MET = pfMET = - all candidates
    // MET (1) in JME-13-003
    retVal.mex = leptonsSum_.mex - pfCandSum_.mex;
    retVal.mey = leptonsSum_.mey - pfCandSum_.mey;
    retVal.sumet = pfCandSum_.sumet - leptonsSum_.sumet;
  }
  if (metType == MvaMEtUtilities::kChHS) {
    //MET = - charged HS
    // MET (2) in JME-13-003
    retVal.mex = leptonsChSum_.mex - pfCandChHSSum_.mex;
    retVal.mey = leptonsChSum_.mey - pfCandChHSSum_.mey;
    retVal.sumet = pfCandChHSSum_.sumet - leptonsChSum_.sumet;
  }
  if (metType == MvaMEtUtilities::kHS) {
    //MET = - charged HS - neutral HS in jets
    // MET (3) in JME-13-003
    retVal.mex = leptonsChSum_.mex - (pfCandChHSSum_.mex + neutralJetHSSum_.mex);
    retVal.mey = leptonsChSum_.mey - (pfCandChHSSum_.mey + neutralJetHSSum_.mey);
    retVal.sumet = pfCandChHSSum_.sumet + neutralJetHSSum_.sumet - leptonsChSum_.sumet;
  }
  if (metType == MvaMEtUtilities::kPU) {
    //MET = - charged PU - neutral PU in jets
    //MET = -recoil in that particular case, - sign not useful for the MVA and then discarded
    //motivated as PU IS its own recoil
    // MET (4) in JME-13-003
    retVal.mex = -(pfCandChPUSum_.mex + neutralJetPUSum_.mex);
    retVal.mey = -(pfCandChPUSum_.mey + neutralJetPUSum_.mey);
    retVal.sumet = pfCandChPUSum_.sumet + neutralJetPUSum_.sumet;
  }
  if (metType == MvaMEtUtilities::kHSMinusNeutralPU) {
    //MET = all candidates - charged PU - neutral PU in jets
    // = all charged HS + all neutrals - neutral PU in jets
    // MET (5) in JME-13-003
    retVal.mex = leptonsSum_.mex - (pfCandSum_.mex - pfCandChPUSum_.mex - neutralJetPUSum_.mex);
    retVal.mey = leptonsSum_.mey - (pfCandSum_.mey - pfCandChPUSum_.mey - neutralJetPUSum_.mey);
    retVal.sumet = (pfCandSum_.sumet - pfCandChPUSum_.sumet - neutralJetPUSum_.sumet) - leptonsSum_.sumet;
  }

  finalize(retVal);
  return retVal;
}

void MvaMEtUtilities::computeAllSums(const std::vector<reco::PUSubMETCandInfo>& jets,
                                     const std::vector<reco::PUSubMETCandInfo>& leptons,
                                     const std::vector<reco::PUSubMETCandInfo>& pfCandidates) {
  cleanedJets_ = cleanJets(jets, leptons, ptThreshold_, 0.5);

  leptonsSum_ = computeCandSum(kLeptons, 0., 0, false, false, leptons);
  leptonsChSum_ = computeCandSum(kLeptons, 0., 0, true, false, leptons);
  pfCandSum_ = computeCandSum(kPFCands, dzCut_, 2, false, false, pfCandidates);
  pfCandChHSSum_ = computeCandSum(kPFCands, dzCut_, 0, false, false, pfCandidates);
  pfCandChPUSum_ = computeCandSum(kPFCands, dzCut_, 1, false, false, pfCandidates);
  neutralJetHSSum_ = computeCandSum(kJets, 0., 0, false, true, jets);
  neutralJetPUSum_ = computeCandSum(kJets, 0., 0, false, false, jets);
}

double MvaMEtUtilities::getLeptonsSumMEX() const { return leptonsSum_.mex; }

double MvaMEtUtilities::getLeptonsSumMEY() const { return leptonsSum_.mey; }

double MvaMEtUtilities::getLeptonsChSumMEX() const { return leptonsChSum_.mex; }

double MvaMEtUtilities::getLeptonsChSumMEY() const { return leptonsChSum_.mey; }

const std::vector<reco::PUSubMETCandInfo>& MvaMEtUtilities::getCleanedJets() const { return cleanedJets_; }