METSignificance.cc

Go to the documentation of this file.

 
// -*- C++ -*-
//
// Package:    METAlgorithms
// Class:      METSignificance
//
//
// Original Author:  Nathan Mirman (Cornell University)
//         Created:  Thu May 30 16:39:52 CDT 2013
//
//
 
#include "RecoMET/METAlgorithms/interface/METSignificance.h"
#include <unordered_set>
 
namespace {
  struct ptr_hash {
    std::size_t operator()(const reco::CandidatePtr& k) const {
      if (k.refCore().isTransient())
        return (unsigned long)k.refCore().productPtr() ^ k.key();
      else
        return k.refCore().id().processIndex() ^ k.refCore().id().productIndex() ^ k.key();
    }
  };
}  // namespace
 
metsig::METSignificance::METSignificance(const edm::ParameterSet& iConfig) {
  edm::ParameterSet cfgParams = iConfig.getParameter<edm::ParameterSet>("parameters");
 
  double dRmatch = cfgParams.getParameter<double>("dRMatch");
  dR2match_ = dRmatch * dRmatch;
 
  jetThreshold_ = cfgParams.getParameter<double>("jetThreshold");
  jetEtas_ = cfgParams.getParameter<std::vector<double> >("jeta");
  jetParams_ = cfgParams.getParameter<std::vector<double> >("jpar");
  pjetParams_ = cfgParams.getParameter<std::vector<double> >("pjpar");
}
 
metsig::METSignificance::~METSignificance() {}
 
reco::METCovMatrix metsig::METSignificance::getCovariance(const edm::View<reco::Jet>& jets,
                                                          const std::vector<edm::Handle<reco::CandidateView> >& leptons,
                                                          const edm::Handle<edm::View<reco::Candidate> >& pfCandidatesH,
                                                          double rho,
                                                          JME::JetResolution& resPtObj,
                                                          JME::JetResolution& resPhiObj,
                                                          JME::JetResolutionScaleFactor& resSFObj,
                                                          bool isRealData,
                                                          double& sumPtUnclustered,
                                                          edm::ValueMap<float> const* weights) {
  //pfcandidates
  const edm::View<reco::Candidate>& pfCandidates = *pfCandidatesH;
 
  // metsig covariance
  double cov_xx = 0;
  double cov_xy = 0;
  double cov_yy = 0;
 
  // for lepton and jet subtraction
  std::unordered_set<reco::CandidatePtr, ptr_hash> footprint;
 
  // subtract leptons out of sumPtUnclustered
  for (const auto& lep_i : leptons) {
    for (const auto& lep : lep_i->ptrs()) {
      if (lep->pt() > 10) {
        for (unsigned int n = 0; n < lep->numberOfSourceCandidatePtrs(); n++)
          footprint.insert(lep->sourceCandidatePtr(n));
      }
    }
  }
 
  std::vector<bool> cleanedJets(jets.size(), false);
  std::transform(jets.begin(), jets.end(), cleanedJets.begin(), [this, &leptons](auto const& jet) -> bool {
    return cleanJet(jet, leptons);
  });
  // subtract jets out of sumPtUnclustered
  auto iCleaned = cleanedJets.begin();
  for (const auto& jet : jets) {
    // disambiguate jets and leptons
    if (!(*iCleaned++))
      continue;
    for (unsigned int n = 0; n < jet.numberOfSourceCandidatePtrs(); n++) {
      footprint.insert(jet.sourceCandidatePtr(n));
    }
  }
 
  // calculate sumPtUnclustered
  for (size_t i = 0; i < pfCandidates.size(); ++i) {
    // check if candidate exists in a lepton or jet
    bool cleancand = true;
    if (footprint.find(pfCandidates.ptrAt(i)) == footprint.end()) {
      float weight = (weights != nullptr) ? (*weights)[pfCandidates.ptrAt(i)] : 1.0;
      //dP4 recovery
      for (const auto& it : footprint) {
        if (it.isNonnull() && it.isAvailable() && (reco::deltaR2(*it, pfCandidates[i]) < 0.00000025)) {
          cleancand = false;
          break;
        }
      }
      // if not, add to sumPtUnclustered
      if (cleancand) {
        sumPtUnclustered += pfCandidates[i].pt() * weight;
      }
    }
  }
 
  // add jets to metsig covariance matrix and subtract them from sumPtUnclustered
  iCleaned = cleanedJets.begin();
  for (const auto& jet : jets) {
    // disambiguate jets and leptons
    if (!(*iCleaned++))
      continue;
 
    double jpt = jet.pt();
    double jeta = jet.eta();
    double feta = std::abs(jeta);
    double c = jet.px() / jet.pt();
    double s = jet.py() / jet.pt();
 
    JME::JetParameters parameters;
    parameters.setJetPt(jpt).setJetEta(jeta).setRho(rho);
 
    // jet energy resolutions
    double sigmapt = resPtObj.getResolution(parameters);
    double sigmaphi = resPhiObj.getResolution(parameters);
    // SF not needed since is already embedded in the sigma in the dataGlobalTag
    //      double sigmaSF = isRealData ? resSFObj.getScaleFactor(parameters) : 1.0;
 
    // split into high-pt and low-pt sector
    if (jpt > jetThreshold_) {
      // high-pt jets enter into the covariance matrix via JER
 
      double scale = 0;
      if (feta < jetEtas_[0])
        scale = jetParams_[0];
      else if (feta < jetEtas_[1])
        scale = jetParams_[1];
      else if (feta < jetEtas_[2])
        scale = jetParams_[2];
      else if (feta < jetEtas_[3])
        scale = jetParams_[3];
      else
        scale = jetParams_[4];
 
      //         double dpt = sigmaSF*scale*jpt*sigmapt;
      double dpt = scale * jpt * sigmapt;
      double dph = jpt * sigmaphi;
 
      cov_xx += dpt * dpt * c * c + dph * dph * s * s;
      cov_xy += (dpt * dpt - dph * dph) * c * s;
      cov_yy += dph * dph * c * c + dpt * dpt * s * s;
 
    } else {
      // add the (corrected) jet to the sumPtUnclustered
      sumPtUnclustered += jpt;
    }
  }
 
  //protection against unphysical events
  if (sumPtUnclustered < 0)
    sumPtUnclustered = 0;
 
  // add pseudo-jet to metsig covariance matrix
  cov_xx += pjetParams_[0] * pjetParams_[0] + pjetParams_[1] * pjetParams_[1] * sumPtUnclustered;
  cov_yy += pjetParams_[0] * pjetParams_[0] + pjetParams_[1] * pjetParams_[1] * sumPtUnclustered;
 
  reco::METCovMatrix cov;
  cov(0, 0) = cov_xx;
  cov(1, 0) = cov_xy;
  cov(0, 1) = cov_xy;
  cov(1, 1) = cov_yy;
 
  return cov;
}
 
double metsig::METSignificance::getSignificance(const reco::METCovMatrix& cov, const reco::MET& met) {
  // covariance matrix determinant
  double det = cov(0, 0) * cov(1, 1) - cov(0, 1) * cov(1, 0);
 
  // invert matrix
  double ncov_xx = cov(1, 1) / det;
  double ncov_xy = -cov(0, 1) / det;
  double ncov_yy = cov(0, 0) / det;
 
  // product of met and inverse of covariance
  double sig = met.px() * met.px() * ncov_xx + 2 * met.px() * met.py() * ncov_xy + met.py() * met.py() * ncov_yy;
 
  return sig;
}
 
bool metsig::METSignificance::cleanJet(const reco::Jet& jet,
                                       const std::vector<edm::Handle<reco::CandidateView> >& leptons) {
  for (const auto& lep_i : leptons) {
    for (const auto& lep : *lep_i) {
      if (reco::deltaR2(lep, jet) < dR2match_) {
        return false;
      }
    }
  }
  return true;
}