NoPileUpMEtUtilities.cc

Go to the documentation of this file.

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

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

#include <algorithm>
#include <cmath>

NoPileUpMEtUtilities::NoPileUpMEtUtilities() {
  minPtDef_ = -1;
  maxPtDef_ = 1000000;
}

NoPileUpMEtUtilities::~NoPileUpMEtUtilities() {}

// namespace NoPileUpMEtUtilities
// {
//-------------------------------------------------------------------------------
// general auxiliary functions
void NoPileUpMEtUtilities::finalizeMEtData(CommonMETData& metData) {
  metData.met = sqrt(metData.mex * metData.mex + metData.mey * metData.mey);
  metData.mez = 0.;
  metData.phi = atan2(metData.mey, metData.mex);
}
//-------------------------------------------------------------------------------

//-------------------------------------------------------------------------------
// auxiliary functions for jets
reco::PUSubMETCandInfoCollection NoPileUpMEtUtilities::cleanJets(
    const reco::PUSubMETCandInfoCollection& jets,
    const std::vector<reco::Candidate::LorentzVector>& leptons,
    double dRoverlap,
    bool invert) {
  reco::PUSubMETCandInfoCollection retVal;
  for (reco::PUSubMETCandInfoCollection::const_iterator jet = jets.begin(); jet != jets.end(); ++jet) {
    bool isOverlap = false;
    for (std::vector<reco::Candidate::LorentzVector>::const_iterator lepton = leptons.begin(); lepton != leptons.end();
         ++lepton) {
      if (deltaR2(jet->p4(), *lepton) < dRoverlap * dRoverlap) {
        isOverlap = true;
        break;
      }
    }
    if ((!isOverlap && !invert) || (isOverlap && invert))
      retVal.push_back(*jet);
  }
  return retVal;
}

CommonMETData NoPileUpMEtUtilities::computeCandidateSum(const reco::PUSubMETCandInfoCollection& cands,
                                                        bool neutralFracOnly,
                                                        double& sumAbsPx,
                                                        double& sumAbsPy) {
  CommonMETData retVal;
  retVal.mex = 0.;
  retVal.mey = 0.;
  retVal.sumet = 0.;
  double retVal_sumAbsPx = 0.;
  double retVal_sumAbsPy = 0.;
  double pFrac = 1;
  for (reco::PUSubMETCandInfoCollection::const_iterator cand = cands.begin(); cand != cands.end(); ++cand) {
    pFrac = 1;
    if (neutralFracOnly)
      pFrac = (1 - cand->chargedEnFrac());

    retVal.mex += cand->p4().px() * pFrac;
    retVal.mey += cand->p4().py() * pFrac;
    retVal.sumet += cand->p4().pt() * pFrac;
    retVal_sumAbsPx += std::abs(cand->p4().px());
    retVal_sumAbsPy += std::abs(cand->p4().py());
  }
  finalizeMEtData(retVal);
  sumAbsPx = retVal_sumAbsPx;
  sumAbsPy = retVal_sumAbsPy;
  return retVal;
}

//-------------------------------------------------------------------------------

//-------------------------------------------------------------------------------
// auxiliary functions for PFCandidates
reco::PUSubMETCandInfoCollection NoPileUpMEtUtilities::cleanPFCandidates(
    const reco::PUSubMETCandInfoCollection& pfCandidates,
    const std::vector<reco::Candidate::LorentzVector>& leptons,
    double dRoverlap,
    bool invert) {
  // invert: false = PFCandidates are required not to overlap with leptons
  //         true  = PFCandidates are required to overlap with leptons

  reco::PUSubMETCandInfoCollection retVal;
  for (reco::PUSubMETCandInfoCollection::const_iterator pfCandidate = pfCandidates.begin();
       pfCandidate != pfCandidates.end();
       ++pfCandidate) {
    bool isOverlap = false;
    for (std::vector<reco::Candidate::LorentzVector>::const_iterator lepton = leptons.begin(); lepton != leptons.end();
         ++lepton) {
      if (deltaR2(pfCandidate->p4(), *lepton) < dRoverlap * dRoverlap) {
        isOverlap = true;
        break;
      }
    }
    if ((!isOverlap && !invert) || (isOverlap && invert))
      retVal.push_back(*pfCandidate);
  }
  return retVal;
}

reco::PUSubMETCandInfoCollection NoPileUpMEtUtilities::selectCandidates(const reco::PUSubMETCandInfoCollection& cands,
                                                                        double minPt,
                                                                        double maxPt,
                                                                        int type,
                                                                        bool isCharged,
                                                                        int isWithinJet) {
  reco::PUSubMETCandInfoCollection retVal;
  for (reco::PUSubMETCandInfoCollection::const_iterator cand = cands.begin(); cand != cands.end(); ++cand) {
    if (isCharged && cand->charge() == 0)
      continue;
    double jetPt = cand->p4().pt();
    if (jetPt < minPt || jetPt > maxPt)
      continue;
    if (type != reco::PUSubMETCandInfo::kUndefined && cand->type() != type)
      continue;

    //for pf candidates
    if (isWithinJet != NoPileUpMEtUtilities::kAll && (cand->isWithinJet() != isWithinJet))
      continue;

    retVal.push_back(*cand);
  }
  return retVal;
}
//-------------------------------------------------------------------------------

//-------------------------------------------------------------------------------
// auxiliary functions to compute different types of MEt/hadronic recoils
//
// NOTE: all pfCandidates and jets passed as function arguments
//       need to be cleaned wrt. leptons
//

void NoPileUpMEtUtilities::computeAllSums(const reco::PUSubMETCandInfoCollection& jets,
                                          const reco::PUSubMETCandInfoCollection& pfCandidates) {
  reco::PUSubMETCandInfoCollection pfcsCh = selectCandidates(
      pfCandidates, minPtDef_, maxPtDef_, reco::PUSubMETCandInfo::kUndefined, false, NoPileUpMEtUtilities::kAll);

  reco::PUSubMETCandInfoCollection pfcsChHS = selectCandidates(
      pfCandidates, minPtDef_, maxPtDef_, reco::PUSubMETCandInfo::kChHS, true, NoPileUpMEtUtilities::kAll);

  reco::PUSubMETCandInfoCollection pfcsChPU = selectCandidates(
      pfCandidates, minPtDef_, maxPtDef_, reco::PUSubMETCandInfo::kChPU, true, NoPileUpMEtUtilities::kAll);

  reco::PUSubMETCandInfoCollection pfcsNUnclustered = selectCandidates(
      pfCandidates, minPtDef_, maxPtDef_, reco::PUSubMETCandInfo::kNeutral, false, NoPileUpMEtUtilities::kOutsideJet);

  reco::PUSubMETCandInfoCollection jetsHS =
      selectCandidates(jets, 10.0, maxPtDef_, reco::PUSubMETCandInfo::kHS, false, NoPileUpMEtUtilities::kAll);
  reco::PUSubMETCandInfoCollection jetsPU =
      selectCandidates(jets, 10.0, maxPtDef_, reco::PUSubMETCandInfo::kPU, false, NoPileUpMEtUtilities::kAll);

  //not used so far
  //_chPfcSum = computeCandidateSum(pfcsCh, false, &_chPfcSumAbsPx, &_chPfcSumAbsPy);
  chHSPfcSum_ = computeCandidateSum(pfcsChHS, false, chHSPfcSumAbsPx_, chHSPfcSumAbsPy_);
  chPUPfcSum_ = computeCandidateSum(pfcsChPU, false, chPUPfcSumAbsPx_, chPUPfcSumAbsPy_);
  nUncPfcSum_ = computeCandidateSum(pfcsNUnclustered, false, nUncPfcSumAbsPx_, nUncPfcSumAbsPy_);

  nHSJetSum_ = computeCandidateSum(jetsHS, true, nHSJetSumAbsPx_, nHSJetSumAbsPy_);
  nPUJetSum_ = computeCandidateSum(jetsPU, true, nPUJetSumAbsPx_, nPUJetSumAbsPy_);
}

CommonMETData NoPileUpMEtUtilities::computeRecoil(int metType, double& sumAbsPx, double& sumAbsPy) {
  CommonMETData retVal;
  double retSumAbsPx = 0.;
  double retSumAbsPy = 0.;

  //never used....
  // if(metType==NoPileUpMEtUtilities::kChMET) {
  //   CommonMETData chPfcSum = computeCandidateSum(_pfcsCh, false, &trackSumAbsPx, &trackSumAbsPy);
  //   retVal.mex   = -chPfcSum.mex;
  //   retVal.mey   = -chPfcSum.mey;
  //   retVal.sumet = chPfcSum.sumet;
  //   retSumAbsPx = ;
  //   retSumAbsPy = ;
  // }
  if (metType == NoPileUpMEtUtilities::kChHSMET) {
    retVal.mex = chHSPfcSum_.mex;
    retVal.mey = chHSPfcSum_.mey;
    retVal.sumet = chHSPfcSum_.sumet;
    retSumAbsPx = chHSPfcSumAbsPx_;
    retSumAbsPy = chHSPfcSumAbsPy_;
  }
  if (metType == NoPileUpMEtUtilities::kChPUMET) {
    retVal.mex = chPUPfcSum_.mex;
    retVal.mey = chPUPfcSum_.mey;
    retVal.sumet = chPUPfcSum_.sumet;
    retSumAbsPx = chPUPfcSumAbsPx_;
    retSumAbsPy = chPUPfcSumAbsPy_;
  }
  if (metType == NoPileUpMEtUtilities::kNeutralUncMET) {
    retVal.mex = nUncPfcSum_.mex;
    retVal.mey = nUncPfcSum_.mey;
    retVal.sumet = nUncPfcSum_.sumet;
    retSumAbsPx = nUncPfcSumAbsPx_;
    retSumAbsPy = nUncPfcSumAbsPy_;
  }
  if (metType == NoPileUpMEtUtilities::kHadronicHSMET) {
    retVal.mex = chHSPfcSum_.mex + nHSJetSum_.mex;
    retVal.mey = chHSPfcSum_.mey + nHSJetSum_.mey;
    retVal.sumet = chHSPfcSum_.sumet + nHSJetSum_.sumet;
    retSumAbsPx = chHSPfcSumAbsPx_ + nHSJetSumAbsPx_;
    retSumAbsPy = chHSPfcSumAbsPy_ + nHSJetSumAbsPy_;
  }
  if (metType == NoPileUpMEtUtilities::kHadronicPUMET) {
    retVal.mex = chPUPfcSum_.mex + nHSJetSum_.mex;
    retVal.mey = chPUPfcSum_.mey + nHSJetSum_.mey;
    retVal.sumet = chPUPfcSum_.sumet + nHSJetSum_.sumet;
    retSumAbsPx = chPUPfcSumAbsPx_ + nHSJetSumAbsPx_;
    retSumAbsPy = chPUPfcSumAbsPy_ + nHSJetSumAbsPy_;
  }

  sumAbsPx = retSumAbsPx;
  sumAbsPy = retSumAbsPy;

  return retVal;
}