text: 
  // -*- C++ -*-
  //
  // Package:    METAlgorithms
  // Class:      METSignificance
  // 
  /**\class METSignificance METSignificance.cc RecoMET/METAlgorithms/src/METSignificance.cc
  Description: [one line class summary]
  Implementation:
  [Notes on implementation]
  */
  //
  // 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");
    
  }
  
  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) {
  
    //pfcandidates
    const edm::View<reco::Candidate>* pfCandidates=pfCandidatesH.product();
  
     // 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 sumPt
     for ( const auto& lep_i : leptons ) {
       for( const auto& lep : *lep_i ) {
         if( lep.pt() > 10 ){
  	 for( unsigned int n=0; n < lep.numberOfSourceCandidatePtrs(); n++ ){
  	   if( lep.sourceCandidatePtr(n).isNonnull() and lep.sourceCandidatePtr(n).isAvailable() ){
  	     footprint.insert(lep.sourceCandidatePtr(n));
  	   } 
  	 }
         }
       }
     }
     // subtract jets out of sumPt
     for(const auto& jet : jets) {
  
       // disambiguate jets and leptons
       if(!cleanJet(jet, leptons) ) continue;
       for( unsigned int n=0; n < jet.numberOfSourceCandidatePtrs(); n++){
         if( jet.sourceCandidatePtr(n).isNonnull() and jet.sourceCandidatePtr(n).isAvailable() ){
  
  	 footprint.insert(jet.sourceCandidatePtr(n));
         }
       }
  
     }
  
     // calculate sumPt
     double sumPt = 0;
     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) {
  	 if( (it->p4()-(*pfCandidates)[i].p4()).Et2()<0.000025 ){
  	   cleancand = false;
  	   break;
  	 }
         }
         // if not, add to sumPt
         if( cleancand ){
  	 sumPt += (*pfCandidates)[i].pt();
         }
       }
     }
     
     // add jets to metsig covariance matrix and subtract them from sumPt
     for(const auto& jet : jets) {
       
       // disambiguate jets and leptons
       if(!cleanJet(jet, leptons) ) 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 sumPt
           sumPt += jpt;
  
        }
  
     }
  
     //protection against unphysical events
     if(sumPt<0) sumPt=0;
   
     // add pseudo-jet to metsig covariance matrix
     cov_xx += pjetParams_[0]*pjetParams_[0] + pjetParams_[1]*pjetParams_[1]*sumPt;
     cov_yy += pjetParams_[0]*pjetParams_[0] + pjetParams_[1]*pjetParams_[1]*sumPt;
  
     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;
  }
  