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();
        }
    };
}

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");
  useDeltaRforFootprint_ = cfgParams.getParameter<bool>("useDeltaRforFootprint");
  
}

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) {
  //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()) {

       //dP4 recovery
       for( const auto& it : footprint) {
     // Special treatment for PUPPI with dR, since jet candidates (puppi) and MET candidates (puppiForMet)
     // can't be matched through the sourceCandidatePtrs and may have different energy, but same direction.
     if((it.isNonnull()) && (it.isAvailable()) &&
            (((!useDeltaRforFootprint_) && ((it->p4()-pfCandidates[i].p4()).Et2()<0.000025)) ||
         (( useDeltaRforFootprint_) && (reco::deltaR2(*it, pfCandidates[i])<0.00000025)))){
       cleancand = false;
       break;
     }
       }
       // if not, add to sumPtUnclustered
       if( cleancand ){
     sumPtUnclustered += pfCandidates[i].pt();
       }
     }
   }

   // 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;
}